{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lead-Acid Models"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare a standard porous-electrode model for lead-acid batteries with two asymptotic reductions. For a more in-depth introduction to PyBaMM models, see the [SPM notebook](./SPM.ipynb). Further details on the models can be found in [[4]](#References)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[33mWARNING: You are using pip version 22.0.4; however, version 22.3.1 is available.\n",
      "You should consider upgrading via the '/home/mrobins/git/PyBaMM/env/bin/python -m pip install --upgrade pip' command.\u001b[0m\u001b[33m\n",
      "\u001b[0mNote: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import pybamm\n",
    "\n",
    "os.chdir(pybamm.__path__[0] + \"/..\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \"Full\" model\n",
    "\n",
    "### Electrolyte Concentration\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial }{\\partial t}\\left(\\epsilon c\\right) = -\\frac{\\partial N}{\\partial x} + sj, \\\\ \n",
    "    N = -\\frac{\\epsilon^b  D(c)}{\\mathcal{C}_\\text{e}}  \\frac{\\partial c}{\\partial x}\\\\\n",
    "    N\\big|_{x=0}= N\\big|_{x=1}=0, \\\\ \n",
    "    c\\big|_{t=0} = 1\n",
    "$$ \n",
    "\n",
    "### Porosity\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial \\epsilon}{\\partial t} = -\\beta^\\text{surf}j, \\\\  \n",
    "    \\epsilon\\big|_{t=0} = \\epsilon^0\n",
    "$$ \n",
    "\n",
    "### Electrolyte Current \n",
    "\n",
    "$$\n",
    "    \\frac{\\partial i_{\\text{e}}}{\\partial x} = j, \\\\\n",
    "    \\mathcal{C}_\\text{e} i_{\\text{e}} = \\epsilon_k^b \\kappa(c) \n",
    "        \\left( \\chi \\frac{\\partial}{\\partial x}\\log(c) \n",
    "        - \\frac{\\partial\\phi_{\\text{e}}}{\\partial x}\\right)\\\\\n",
    "    i_{\\text{e}}\\big|_{x=0}= i_{\\text{e}}\\big|_{x=1}=0,\n",
    "$$\n",
    "\n",
    "### Electrode Current \n",
    "\n",
    "$$\n",
    "    \\frac{\\partial i_{\\text{s}}}{\\partial x} = -j,\\\\\n",
    "    i_{\\text{s}} = -\\sigma\\frac{\\partial\\phi_{\\text{s}}}{\\partial x},\\\\\n",
    "    \\phi_{\\text{s}}\\big|_{x=0}\n",
    "        = i_{\\text{s}}\\big|_{x=l_\\text{n}} \n",
    "        = i_{\\text{s}}\\big|_{x=1-l_\\text{p}} = 0, \\\\\n",
    "    i_{\\text{s}}\\big|_{x=1}=\\mathcal{I},\n",
    "$$\n",
    "\n",
    "\n",
    "### Interfacial current density\n",
    "\n",
    "$$\n",
    "j = \\begin{cases}\n",
    "    2j_0(c) \\sinh\\left(\\eta\\right), \\quad &0 < x < l_\\text{n} \\\\\n",
    "    0, \\quad &l_\\text{n} < x < 1-l_\\text{p} \\\\\n",
    "    2j_0(c) \\sinh\\left(\\eta\\right), \\quad &1-l_\\text{p} < x < 1 \n",
    "\\end{cases}\n",
    "\\\\\n",
    "\\eta = \\phi_{\\text{s}} - \\phi_\\text{e} - U(c),\n",
    "$$\n",
    "\n",
    "\n",
    "This model is implemented in PyBaMM as the `Full` model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "full = pybamm.lead_acid.Full()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \"Leading-order\" model \n",
    "\n",
    "$$\n",
    "    \\frac{\\mathrm{d} }{\\mathrm{d} t}\\left(\\epsilon c\\right) = (s_\\text{n} - s_\\text{p})\\mathrm{I}, \\\\\n",
    "    \\frac{\\mathrm{d} \\epsilon}{\\mathrm{d} t} = -\\beta^\\text{surf}j, \\\\\n",
    "    j = \\begin{cases}\n",
    "        \\mathrm{I}/l_\\text{n}, \\quad &0 < x < l_\\text{n} \\\\\n",
    "        0, \\quad &l_\\text{n} < x < 1-l_\\text{p} \\\\\n",
    "        -\\mathrm{I}/l_\\text{p}, \\quad &1-l_\\text{p} < x < 1 \n",
    "    \\end{cases} \\\\\n",
    "    \\phi_\\text{e} = -U_\\text{n}(c) + \\sinh^{-1}\\left(\\frac{\\mathrm{I}}{2l_\\text{n}j_{0\\text{n}}(c)}\\right) \\\\\n",
    "    V = -\\phi_\\text{e} + U_\\text{p}(c) - \\sinh^{-1}\\left(\\frac{\\mathrm{I}}{2l_\\text{p}j_{0\\text{p}}(c)}\\right) \\\\\n",
    "$$\n",
    "\n",
    "This model is implemented in PyBaMM as `LOQS` (leading-order quasi-static)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "loqs = pybamm.lead_acid.LOQS()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the models"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We load process parameters for each model, using the same set of (default) parameters for all. In anticipation of changing the current later, we make current an input parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load models\n",
    "models = [loqs, full]\n",
    "\n",
    "# process parameters\n",
    "param = models[0].default_parameter_values\n",
    "param[\"Current function [A]\"] = \"[input]\"\n",
    "for model in models:\n",
    "    param.process_model(model)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we discretise the models, using the default settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "for model in models:\n",
    "    # load and process default geometry\n",
    "    geometry = model.default_geometry\n",
    "    param.process_geometry(geometry)\n",
    "\n",
    "    # discretise using default settings\n",
    "    mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)\n",
    "    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)\n",
    "    disc.process_model(model)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we solve each model using CasADi's solver and a current of 1A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solved the LOQS model in 64.611 ms\n",
      "Solved the Full model in 736.504 ms\n"
     ]
    }
   ],
   "source": [
    "timer = pybamm.Timer()\n",
    "solutions = {}\n",
    "t_eval = np.linspace(0, 3600 * 17, 100)  # time in seconds\n",
    "for model in models:\n",
    "    solver = pybamm.CasadiSolver()\n",
    "    timer.reset()\n",
    "    solution = solver.solve(model, t_eval, inputs={\"Current function [A]\": 1})\n",
    "    print(f\"Solved the {model.name} in {timer.time()}\")\n",
    "    solutions[model] = solution"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To plot the results, the variables are extracted from the solutions dictionary. For example, we can compare the voltages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAG2CAYAAACNhdkhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB8iElEQVR4nO3dd3QUVRvH8e9uOiEJPaEm9A5BehGQXqR3qSJNAUVFARtYaSqiIIj0JohKR5Aaeu9VQDqEThJIT+b9Y2V5IwmBkGRTfp9z9hzmzp2ZZ4YtT+7cuddkGIaBiIiIiMTJbOsARERERFI6JUwiIiIi8VDCJCIiIhIPJUwiIiIi8VDCJCIiIhIPJUwiIiIi8VDCJCIiIhIPe1sHkJpER0dz9epV3NzcMJlMtg5HREREnoJhGAQFBZErVy7M5oS1FSlhegZXr14lb968tg5DREREEuDSpUvkyZMnQdsqYXoGbm5ugOWCu7u72zgaEREReRqBgYHkzZvX+jueEEqYnsHD23Du7u5KmERERFKZ5+lOo07fIiIiIvFQwiQiIiISDyVMIiIiIvFQwiQiIiISDyVMIiIiIvFQwiQiIiISDyVMIiIiIvHQOEwiIilYREQEUVFRtg5DJEWxs7PDwcEhWY+phElEJAUKDAzk1q1bhIWF2ToUkRTJycmJbNmyJdtA0kqYRERSmMDAQK5cuULGjBnJli0bDg4OmvBb5F+GYRAREUFAQABXrlwBSJakSQmTiEgKc+vWLTJmzEiePHmUKInEwsXFBTc3Ny5fvsytW7eSJWFSp28RkRQkIiKCsLAwPDw8lCyJPIHJZMLDw4OwsDAiIiKS/HhKmFKAawEhtJm0nX0X7tg6FBGxsYcdvJO7Q6tIavTwc5IcD0YoYUoBRi3Zy8CrQ/l08lw+XHyEgJCkz5RFJGVT65JI/JLzc6KEycaCgkPoenk4te0OscDxC67tWULdb/xYfugqhmHYOjwRERFBCZPNuZkjKJ/TCYAMpjB+dviGhiErGfjLAXrM2MPF28E2jlBERESUMNmaszumrouhZCsA7EwGXzpM5337BWz++zoNvvNj0qazRERF2zhQERGR9EsJU0rg4AxtpkO1N61Fb9gvY5zDj0RHhDF69Ule/n4r+y7ctWGQIiK2YzKZnrm/yuXLlxk8eDClS5fG3d0dV1dXihQpwuuvv87x48fj3T4yMpIpU6ZQr149PD09cXJyImfOnDRr1oyFCxc+sdtEWFgY48aNo2rVqnh4eODo6EjOnDmpUKECAwcOZOXKlc90LimdyWTCx8fnufczc+ZMTCYTI0aMeO59JTYlTCmF2QwNPocmXwOWL4WWdttZ4fgBvqYznLoeRNvJ29UpXETkKcybN48iRYrwzTffcPfuXerWrUuTJk0AmDx5MmXKlGHMmDFxbn/x4kXKlStH37592bp1K6VLl6ZNmzYULlyY1atX07FjR+rUqcO9e/ce2/bevXtUq1aNd955h3379uHr60ubNm2oUKECV69eZcKECbz77rtJdeqSRDRwZUpTqTe45YTfX4PIUIqYr/C70wimRTbm28i2zNt1kTXHrjO8WQleLpNTT9KIiPzHkiVL6Nq1Kw4ODvz000/07t07xnflsmXL6NatG0OGDMHFxYWBAwfG2D4gIIBatWpx/vx5WrduzU8//US2bNms6y9evEjXrl3ZtGkTTZo0YfPmzdjbP/o5/eSTT9i/fz++vr6sWLGC3Llzx9j/7t27Wb16dRKdvSQVtTClRMVfht4bIKcvAHZE08d+JWuchlHJdIJb98OsncIv3VGncBGRh+7fv0/v3r0xDIMpU6bQp0+fx/6wbN68OYsXL8ZkMvHee+9x4cKFGOuHDBnC+fPnqVu3Lr/++muMZAkgX758/Pnnn5QoUYIdO3Ywbty4GOt///13AMaOHftYsgRQqVIlPvnkk8Q4XUlGSphSKs+S0Gs91BsBdpan6LxN/vzq9Dlf2f+MO/fx+/sm9cepU7iIyEOzZ8/m1q1bVKxYke7du8dZ76WXXqJ169aEhYUxceJEa/nt27eZNWsWAOPHj8fOzi7W7TNkyMDo0aMB+O6772IMnHjz5k0AsmfP/tznA9CjRw9MJhObNm1i3bp11KxZEzc3N3LkyEHv3r0JCAgA4MaNG/Tt25fcuXPj7OxMpUqV2LRpU5z7nTNnDjVq1MDd3Z0MGTJQpkwZRo4cSWhoaKz179y5w4ABA8iVKxfOzs6UKFGC8ePHxzsEzq5du2jXrh05c+bE0dGRPHny0KtXLy5evJjga2ILSphSMjt7qPE29NsKeStbi1+x38hG5/d52byD0IgoRq8+SbMf1ClcRORhZ+pXXnkl3rqdO3cG4M8//7SWbdiwgdDQUMqWLUvJkiWfuH3jxo3JnDkzV69e5fDhw9byvHnzApa+Uok5nt7ixYtp1KgRhmHQqFEjnJycmDp1Ki1atODWrVtUrVqVNWvW8OKLL+Lr68uePXto1KgRR44ceWxfffv2pVu3buzbt48XX3yRpk2bcu3aNT744APq1KlDcHDMuxd3796lRo0aTJw4EcMwaNGiBblz52bw4MG8+eabj+3/oR9//JFq1arxxx9/4O3tTcuWLcmaNSvTpk2jQoUKnDhxItGuT5Iz5KkFBAQYgBEQEJD8B4+KNIydPxnGl7kNY7i79bXho5pGtSEzDO8hKwyfoSuMDxcfNu4Fhyd/fCKSKEJCQozjx48bISEhtg4lRQGMp/nJyp07twEYW7ZsibfuhQsXDMAwmUxGeLjle/PDDz80AOO11157qrheeuklAzCmTZtmLRs5cqQ13mLFihlDhw41Fi9ebFy6dOmp9vlf3bt3NwDDbDYbK1assJYHBgYapUqVMgCjRIkSRpcuXaznYRiG8dFHHxmA0a1btxj7++233wzAyJUrl/H3339by+/du2fUqFHDAIx33303xjb9+vUzAKNRo0bGgwcPrOW7du0yMmbMaACGt7d3jG127Nhh2NnZGblz5zb27t0bY93UqVMNwKhcuXKM8hkzZhiAMXz48Ke6Nk/7eUmM3291+k4tzHZQuQ8Uawqr3oNTlr+iXrI7yHq79/k6oi0zohoxd+ejTuFNS6tTuEha1OyHrdwMCrN1GE8lu5sTywfWSLbj3b5923Lcp7gdliNHDgAMw+DOnTt4eno+0/b/v49bt25Zy95//33u3r3LuHHjOHnyJKNGjbKuK1myJAMGDKBPnz6Yzc92k+eVV16hadOm1mU3Nzd69+7NW2+9xeXLl9m6dWuMOQgHDx7Ml19+iZ+fX4z9fP/99wAMHz6cwoULW8s9PDyYOHEivr6+/PTTT3zxxRc4Ozvz4MEDZs2ahdlsZsKECWTIkMG6TaVKlejfv7/19uT/GzVqFFFRUUyePJny5cvHWPfaa6+xbNkyli1bxoEDByhXrtwzXQtbUMKU2njkhk7z4cRyS+IUdA1nwvjIYR6t7bcxJLwXR4IKMGD+AX4vepnPWpQib5YM8e9XRFKNm0Fh+AfG3s9EbOP/+zCZzWZGjx7NW2+9xW+//cbmzZvZs2cPFy9e5NixY7z++uusWbOG33///ZmSpgYNGjxWVqBAAQAqVKhA5syZY6zz8PAgS5YsXLt2zVoWERHBzp07gUe3JP9fmTJlKFOmDIcOHeLgwYNUqVKFffv2ERISQqVKlShYsOBj23Tq1OmxhCk6Opr169eTIUMGGjZsGOv5vPjiiyxbtozdu3crYZIkVLwZ5K8FGz6H3T8DBiVM51nq9AmzIuvzdWR7Np66SYNxm3m7fmFerZ4fBzt1WRNJC7K7Odk6hKeW3LFmzZqVK1eucPPmTYoWLfrEujdu3AAsgy5myZLFuj086rgdn4f7+O+TdAC5cuXizTfftPbxOXHiBF9//TXTp09nyZIl/PLLL7EmLXGJ7Ym7jBkzxrnu4fqHrWZgaYELDw8nW7ZsuLq6xrqNj48Phw4d4sqVKwBcvXoVAG9v7zjr/9etW7e4f/8+AI6OjnGc0aO6qUGKTJhGjhzJH3/8wcmTJ3FxcaFatWqMHj36iW/+Y8eO8cknn7Bv3z4uXLjAuHHjGDRoUIw6I0aM4NNPP41RVrRoUU6ePJkUp5H0nN2hyVgo0wGWvwXXj2Immlft19DEfi8fhXdnbUQFvlp1ksUHrjKydWl882ayddQi8pyS8xZXalO2bFmuXLnC3r17qVHjyddp//79gOU22cNbWWXLlgVg79698R4rOjqaQ4cOWfcRn+LFizNt2jTu3r3L4sWLWbly5TMlTE9qjXrW23tPkhhdOaKjLU9uZ8yYkTZt2jyx7tNcu5QgRSZMfn5+9O/fn4oVKxIZGckHH3xAgwYNOH78eJwZcXBwMAUKFKBdu3a8/fbbce67ZMmSrFu3zrr8/4ONpVp5KkCfTbBjImwaBZEheHKbnx2/ZU1UBYZHdOfENWj14za6V/Xh3QZFcHN2iHe3IiKpTZMmTVi1ahW//PLLY380/9f8+fMBy9NuD9WpUwcnJycOHTrE8ePHKVGiRJzbr169mjt37pA1a1YqVar01DHWqVOHxYsX26RlJWvWrDg6OnLr1i0ePHgQ62/q+fPngUetVjlz5gR4bLyqh2Irz5YtG87OzpjNZmbMmJEm+tOmyHs0q1evpkePHpQsWZKyZcsyc+ZMLl68yL59++LcpmLFiowdO5aOHTvi5BR3E7C9vT1eXl7WV2zNqKmSnQPUGAT9d0LButbihnZ72eD8Pt3t1mAyopm5/Tz1v93M6qP+totVRCSJdO/enaxZs7J7927mzJkTZz0/Pz9+++03nJyc6N+/v7U8a9as1vGbBg0aZG0p+a+QkBDef/99AN58880Yf3wb8QwlcObMGSDu22hJycHBgSpVqgCwYMGCx9YfPXqUQ4cOkTFjRnx9fQEoX748Li4u7Nu3j3/++eexbWLbj729PbVr1yYwMJD169cn7knYSIpMmP7r4aBcD+8xP4/Tp0+TK1cuChQoQOfOnZ84cFZYWBiBgYExXileZh/o8ju0mQaulqc8MhDCpw6zWOw0nBKm8/gHhtJv7j56z97L1Xshto1XRCQRZcyYkalTp2IymejVqxfTpk17LIFZsWIFLVu2xDAMxowZ81jfnNGjR+Pj48PatWvp0KFDjD5AAJcuXaJJkyYcO3aM6tWr88EHH8RYX61aNWbMmMGDBw8ei2/FihVMnjwZgLZt2ybGKT+zh1PBjBgxIkYCFBQUxIABAzAMg759++Ls7AxYrmnXrl2Jiopi4MCBhIQ8+t3Yu3cvEyZMiPU4H374IWazmVdffTXWATTv37/P9OnTY+wvRUvwgATJJCoqymjatKlRvXr1p97G29vbGDdu3GPlq1atMn799Vfj0KFDxurVq42qVasa+fLlMwIDA2Pdz/Dhw61jafz/yybjMCVE8B3DWDowxrhNkcMzGZM/7GQUG/Kb4T1khVHi4z+NaVv+MSKjom0drYgYGocpLg+/fytXrhzn6+eff7bWnzNnjuHi4mIARp48eYzWrVsb7dq1M4oUKWIAhp2dnTFq1Kg4j3fhwgWjZMmSBmA4Ozsb9erVMzp16mTUqlXLsLe3NwCjTp06xrVr1x7b1sPDw7pdtWrVjI4dOxotWrQwihUrZj2Pfv36PfW5PxyHaePGjY+t27hxowEY3bt3j3Vbb2/vWMev6tOnjwEYLi4uRtOmTY127doZ2bNnNwCjSpUqMcZaMgzDuH37tlG0aFEDMHLmzGl06NDBaNCggWFvb2/0798/1nGYDMMwJk2aZNjZ2RmAUapUKaN169ZGhw4djMqVKxtOTk4GYNy9e9daPyWPw5TiE6Z+/foZ3t7ezzTgV1wJ03/dvXvXcHd3N6ZOnRrr+tDQUCMgIMD6unTpUupKmB46v90wfqgYI3G6Mryg0X3Y54b3kBWG95AVxsvfbzGOXL5n60hF0j0lTLGL7Y/X/77++yN78eJF49133zVKlChhuLq6Gi4uLkbBggWNvn37GkePHo33mOHh4cZPP/1k1KlTx8iWLZthNputx2rbtq0RHR37H5qHDh0yRo8ebTRo0MAoWLCgkSFDBsPJycnIly+f0bZtW2PVqlXPdO5JkTAZhmHMnj3bqFatmpExY0bD2dnZKFmypPHll18awcHBsda/deuW8frrrxteXl6Gk5OTUaxYMePrr782oqOj40yYDMMwDhw4YHTv3t3w9vY2HB0djUyZMhklS5Y0evbsaaxYsSLGdUzJCZPJMBJx3PZENmDAAJYuXcrmzZvJnz//U2/n4+PDoEGD4u3wB5a+T/Xq1WPkyJHx1g0MDMTDw4OAgADc3d2fOp4UITIcto2HzWMh6tGAd8ujqvBZRDdukgmzCXpWz8/b9Yvg6pQGOsOLpEKhoaGcO3eO/PnzW2+JSMqxYcMGGjduTFRUFIsXL6ZZs2a2Dilde9rPS2L8fqfIPkyGYTBgwAAWL17Mhg0bnilZehb379/n7Nmz1icA0jR7R6j1Hry+HXxetBY3s9vJBuf36GS3HsOIZurWczQYt5kNJ6/bMFgRkZSpTp06zJs3D8Mw6NChAzt27LB1SJJMUmTC1L9/f+bOncv8+fNxc3PD398ff3//GB3DunXrxrBhw6zL4eHhHDx4kIMHDxIeHs6VK1c4ePCg9WkEsAwT7+fnx/nz59m+fTutWrXCzs6OTp06Jev52VS2QtB9ObScBC6WTvRuPGCkwzQWOX1OIdNlrtwLoefMvfSft58bGk1YRCSGtm3b8ttvv/H++++zd+/eRJ1gV1KuFHlLLq7xGmbMmEGPHj0AqF27Nj4+PsycOROwjBsRW0tUrVq1rL3zO3bsyObNm7l9+zbZs2enRo0afPnll7EO9R6bVH1LLjYPbsNfH8KhX6xFkdjzY+TLTIxsSRiOuDnbM6RRMV6plA+zOfWPoyGS0umWnMjTS85bcikyYUqp0lzC9NA/m2DF23Dn0eOlF8jJ0PCe7Ii2jMBa3jszI1uXpoinm42CFEkflDCJPL1034dJklmB2pa+TS++C2ZLZ29vrvGL45eMsf+JTASx78Jdmozfwtg1JwmNiHry/kRERNIYJUxi4eACdT+Bvlsgz6Mh/tvb+7HJ+T1amrcSGR3NxI1nafTdZrafSR2TJYqIiCQGJUwSk2cJ6LkGmn4LTpZmy0wE8p3jj8x1HEU+03XO3w7mlam7eOfXg9x5EG7jgEVERJKeEiZ5nNkMFV+D/ruhRAtrcQ3zEdY6DaGf3TLsieSP/Veo+80mft93WU+JiIhImqaESeLmnhPaz4ZOC8A9DwBOhDPUYQErnT/G13SGu8ERvLvoEF2m7eL8rcfnTRIREUkLlDBJ/Io2hv47ocobYLK8ZYpygT+chjPCfiYZCWbbmds0/G4zEzeeISIq9tm9RUREUislTPJ0nNyg0UjotR68SgNgxqCH/V9scB5CA/MewiKjGbvmFC9/v5V9F+7aOGAREZHEo4RJnk3uF6D3Jqj/OThkACAHt5niOI6fHMbhyR1OXQ+i7eTtfLzkKIGhEbaNV0REJBEoYZJnZ2cP1d+EN3ZAoXrW4oZ2e9jo8h5d7NaCEc2cnReo/60fq4/62zBYERGR56eESRIusw90/g3aTAPX7ABkMEL4wmEGfzh9RhHTJa4HhtFv7j56z97LtYCQJ+9PRCQOJpPpia/atWs/9/59fHxilJ0/fz5R9p2SJOY5jRgxApPJZJ2iLK2zt3UAksqZTFC6LRSsA2s/gQNzAChn+ptVTh8wKbIZEyJbsvb4dbafucV7DYvStaoPdpqXTkQSoHv37rGWFytWLJkjkfRGCZMkjgxZoMUEKNMBVgyC22ewJ4qB9ktobr/LMi9deElGLD/O4oNXGdmqNCVypaH5+EQkWaSX1gxJeXRLThJX/heh3zao+T6YHYDH56U7dOkezSZsZdSfJwkJ17x0IiKS8ilhksTn4Ax1PoR+WyBvZWtxe3s/Nrq8T3PzNqKio5nsd5aG321my+mbNgxWRNKSmTNnYjKZGDFiRKzra9eujclk4vz580ly/P/f/8KFC6lYsSIZMmQgd+7cvP/++4SHW6aTOnv2LJ06dSJHjhxkyJCBl156icOHD8e6z8jISH744QfKly9PxowZyZgxI5UqVWLSpElERcX+R+elS5fo2rUr2bNnJ0OGDJQvX565c+fGG//q1atp2rQp2bNnx8nJiQIFCvDOO+9w+/bthF+UNEIJkySdHMXh1dUx5qXLbATwveNEZjuOIY/pBhfvBNN12m7eWXiQ2/fDbBywiEjiGD9+PF26dCFTpkw0atSI8PBwxo4dS+/evTl9+jRVqlTh4MGD1KlTh0KFCrFp0yZeeuklrl+/HmM/UVFRtGjRgjfffJMzZ85Qv3596tWrx8mTJ3njjTdo164d0dExBws+d+4clSpVYu7cubi7u9OiRQtcXV3p1q0b33zzTZwxDx06lMaNG7Nu3TqKFi1K8+bNsbe3Z9y4cVSuXPmx2NIbJUyStP5/Xrriza3FNc2HWOc8hF52K7Ejij8OXKHet36al05E0oSpU6eyY8cO1q5dyx9//MGRI0fw9PRkzpw5tGjRgl69enH8+HEWLFjAoUOH6Nq1K3fu3OHHH3+MsZ/vvvuOVatWUbJkSf7++28WL17MkiVLOHXqFEWLFmXx4sWPbfPGG2/g7+9Pz549OXXqFL/88gubN29m6dKlTJo0KdZ4Fy1axOjRoylVqhTHjh1j69atLFq0iFOnTvHJJ59w9uxZ3nrrrSS7XqmBOn1L8nDPCR3mwMmVsHIwBF3F2QjjI4d5tHbYznthvTgWnJ93Fx1i8YErfNmqFN5ZXW0dtUjK9FMtuH/D1lE8nYw5oK9fou3OZIr9Cdtz5849NiyALQ0aNIgKFSpYl728vHjllVcYN24cYWFhfPbZZ9ZzMZlMDB48mDlz5uDnF/Naff/99wB8++23eHp6Wstz5szJ2LFjad68OePHj2fAgAEA/PPPP6xevRp3d3e+/fZb7O0f/cw3a9aMtm3bsnDhwsfi/fLLLwH45ZdfKFSokLX84e3NZcuW8dtvv3Hr1i2yZcv2vJcnVVLCJMmrWFPweRE2fA67fwYMSnCO5U4fMzWyMeMi27D1zC0ajNvMoHpF6PVifhzs1BAqEsP9GxB01dZR2ERcwwpkzJgxmSN5sgYNGjxWVqBAAcDSz8nBwSHWddeuXbOWXbx4kYsXL5I9e/ZY9/fyyy+TKVMmzpw5g7+/P15eXmzduhWARo0a4eHh8dg2nTp1eixhunHjBocOHaJw4cKUKlXqsW1MJhPVq1fn4MGD7Nu3j4YNG8Z3+mmSEiZJfs7u0GQslG4Py9+EG8cxE00f+5W8bL+HoeGvsjmyLKNXn2TZoauMal2asnkz2TpqkZQjYw5bR/D0EjnW1DKsQO7cuR8re5jUPWldWNijvpxXr1qSYm9v71iPYTKZ8Pb25t69e1y5cgUvL694t4mtFe5hB/jTp0/H2YL30K1bt564Pi1TwiS2k7ci9N0M28aD3xiICiMXN5jtOJqlUdX4LKIrJ65Bqx+30b2aD+82KEpGJ71lRRLzFld6898O0knFbI67ZfxJ655VfAnO03h4Tby8vOJtPYorEUsP9OsjtmXnADUHQ8lWsPwtOL8FgBZ223nJ/gifhb/Cb1E1mbHtPGuO+vN5y1LULe4Zz05FJL1ydHQE4P79+7Guv3TpUnKG81xy5coFwIULF+Ks83Ddw1arnDlzPnGb2Mrz5MkDQLZs2VJNC54tqHOIpAxZC0L35dBiIjhnAsDdCOJrh5+Y7zgSb5M/VwNCeW3WXvrP38+NoFDbxisiKdLDhOHvv/9+bN3ff//NxYsXkzukBMuXLx/58uXj5s2brF+//rH1K1eu5O7duxQqVAgvLy8AatSoAVjGUwoMDHxsmwULFjxWlidPHooVK8bx48djvW5ioYRJUg6TCcp1gQF7oVRba3E181HWOg2ln90y7Ilk5eFr1PvGjwW7L2oIAhGJ4eFAkX/++Sf79u2zlt+6dYtevXol2y25xDJw4EAA3nnnHW7efDTIr7+/P++99x5AjMf9CxYsSIMGDQgMDOTdd9+NMbDlqlWrWLRoUazH+fjjj4mOjqZNmzYcPHjwsfW3b9/m559/ToxTSrWUMEnKkzE7tJ0GrywCj7wAOBLOUIcFrHT+mDKmswSGRjL0jyN0nLKTf27G3vQuIulPxowZGTx4MJGRkdSoUYNGjRrRuHFjihQpQlRUFFWrVrV1iM/k7bffpnHjxhw+fJjChQvTunVrWrVqRZEiRThx4gQtW7bkjTfeiLHNpEmT8PT0ZOrUqRQtWpROnTpRq1YtXn75Zfr27RvrcV555RU++OADjh49Svny5Slfvjzt27enXbt2vPDCC3h6eloTtPRKCZOkXEUawBs7oUp/MFneqkW5wBKn4XxsP4cMhLLr3B0ajd/ChA2nCY9MXX85ikjSGDFiBGPHjiVPnjxs2LCBo0eP0rNnT9auXWvt45Ra2NnZsWzZMsaPH0+BAgVYs2YNf/31F0WLFmXixIn89ttvj3UiL1CgALt27eKVV17h3r17LFmyhMDAQGbMmMHgwYPjPNaXX36Jn58fbdq0wd/fnyVLlrBx40aioqJ4/fXXWbZsWVKfbopmMnRP46kFBgbi4eFBQEAA7u7utg4nfbmyD5a9BdePWIuumbIzLOxVNkX7AlDU042RbUrzQr7MNgpS5PmFhoZy7tw58ufPj7Ozs63DEUnRnvbzkhi/32phktQhd3nosxHqfQr2lg9FTuMmMx3H8L3DBLISwKnrQbSZtJ1Plh4lKDTCxgGLiEhaooRJUg87B6gxCN7YAflrWYub221no8v7tDFvxjAMZu+4QINxm1l3PH1PFCkiIolHCZOkPlkKQLel0OLHGEMQfOM4mflOo8hnus61gFB6zd5L/3kagkBERJ6fEiZJnUwmKNf58SEITEdY6zyUvnbLsSOKlUc0BIGIiDw/JUySuv3/EATultFqnYwwhjn8wgrnTyhpOqchCERE5LkpYZK0oUgD6L8TKr8OWOZWKs45ljt9zAf283DREAQiIvIclDBJ2uHkBo1HQa91kKMEAGai6WO/kvUuw6hhPkJ4ZDRf//U3zX7YyoGLd20csEjcdAtZJH7J+TlRwiRpT54K0McP6nwEdk4A5DKuM9dxJF87TCYTQZy6HkTrSdsZsewYD8IibRywyCMPByH8/yktRCR2Dz8n/x28MykoYZK0yd4Rar4Hr28H7+rW4rZ2m9nk8j7NzdsxDIOZ28/TYNxmNp68YcNgRR5xcHDAzs6OkJAQW4cikuKFhIRgZ2eHg4NDkh9LCZOkbdkKQfcV0Gw8OHkAkMkI4HvHCcxw/Jpc3OLKvRBenbmHN385wK37YTYOWNI7k8lEhgwZCAgIUCuTyBNERUUREBBAhgwZMJlMSX48TY3yDDQ1SioX5A+r3oMTj+ZDCjG5MCq8PXOi6hONmUwZHPioaQnavJA7WT6AIrEJDw/n/Pnz2NvbkyVLFpycnPR+FPmXYRiEhYVx584dIiMj8fHxiXeOwMT4/VbC9AyUMKURJ1bAqsEQdM1adIjCDA7rzWnDMjRBjULZ+KpVafJlzWCrKCWdCw4O5tatWzx48MDWoYikSK6urmTLlo0MGeL/nlbClMyUMKUhoQGwbgTsnW4tisSeiZHNmRjZgnAccHYw8079IvSsnh97O929FtuIjIwkMlIPJoj8P3t7e+zt7Z+6vhKmZKaEKQ26sB2WvQm3T1uLzpny8G5oL/YbRQAoldud0W3KUDKXh62iFBGR55AYv9/6s1nSN+9q0G+r5Yk6s+WvlfzGZX53+pTPHWaQkWCOXgmk+YRtjPrzJKER6oQrIpIeqYXpGaiFKY3zPwrLBsLV/daim6ZsDAnrwYboFwDwyZqBr1qXplrBbLaKUkREnlGy3ZKbPXt2gnYem27duiXavpKbEqZ0IDoKdv0EGz6HiGBr8croqnwS3o3bWG7LdayYl2GNi+ORIenH/hARkeeTbAmT2WxOtEdaU/O4IkqY0pG752HF23B2g7UoyOTGiLDO/B79ImAiW0YnPmtRksalvPTIt4hICpasCZOvry8tWrRI0EEAlixZwuHDh5UwSephGHB4IaweCiGP5p3bbpTh/fCeXDZyANCghCeftSiFl4ezrSIVEZEnSIzf76d+Js/X15fhw4cn6CAA58+f5/DhwwneXiTZmUxQtiMUrGtJmo7+BkA102HWOw9lTHhbZkQ14q/j19lx9jZDmxSjU8V8mM1qbRIRSWue6ik5d3f3pxoY6klcXFzUKiOpU8bs0HYavPIruFsGtnQyQvnYYS7LnUdQ1HSRoLBIPlx8lI4/7+Tszfs2DlhERBKbnpJ7BrolJ4QFwbpPYc9UwPLRicKOiZHNmBDZinAccLQ381bdwvSpWQAHDXgpImJzGrgymSlhEquLOy1DENz621p03pSHd/5vwMviOd0Z3aY0ZfJkslGQIiICyThwZbt27Vi+fLmG5xd5KF+Vfwe8fN864KXPvwNefuowE1dCOHEtkJYTt/HlyuOEhKfehx1EROQZhxXIkiULHTp0oEuXLlSpUiU54ktR1MIksbp+zNLadGWfteiGKRvvh/VkU7QvAPmyZGBU69JUK6QBL0VEkluy3ZIbM2YM8+bN48iRI5aNTCYKFixI165d6dy5MwUKFEjQwVMbJUwSp+go2DUZNnwRY8DLpdE1GBHehbtY3i8dKuTlgyYa8FJEJDklex+mI0eOMGfOHBYsWMDly5etg/VVrVqVrl270r59ezJnzpygQFIDJUwSr7vnYflb8M8ma1GgyYOPwrqyLLoqYCK7mxOftyhJo1I5bRWliEi6YtNO35s2bWLOnDn88ccfBAQEYDKZcHBwoEmTJnTt2pWXX34ZB4e09Ve0EiZ5KoYBB+fBmg8gNMBa7Ge8wNCwV7lGVgAalfTisxYlyeGuAS9FRJJSinhKLiwsjOXLlzN37lxWr15NeHg4JpOJTJky0b59ezp37kyNGjWe5xAphhImeSZB12HVYDixzFoUYsrAl+EdmBdVFwMz7s72fNi0OO0r5NX0KiIiSSRFJEz/7+7duyxcuJD58+ezbds2ywFMpjTzdJ0SJkmQE8th5btw/7q1aD/FGRz2Gv8YuQCoVjArI1uXxjurq62iFBFJs5JtWIGnlTlzZqpXr07VqlXJkiULhmGgYZ4k3SveDPrvgnJdrUUvcII1zh/wut0y7Ilk+9nbNPxuMz9v/ofIqGgbBisiIrFJlBamK1euMH/+/BhP0hmGQalSpejatSvvvffecweaEqiFSZ7bP36w/E1L5/B/nTIV4J3QXhwzfAAok8eD0W3KUDyn3mMiIonBprfkgoKCWLRoEXPnzmXz5s3W1qScOXPSqVMnunbtStmyZRMUVEqlhEkSRXgwbPwSdv4IhqU1KRo7JkU25fvI1oThiL3ZxOu1CzKgTiGc7O1sHLCISOqW7AlTZGQkK1euZN68eaxYsYKwsDAMw8DV1ZVWrVrRtWtX6tati9mcNufPUsIkieryPlg2AG4ctxZdNOfm7ZBe7DOKAlAoR0ZGtylNee8stopSRCTVS7aEadu2bcydO5dFixZx9+5dDMPAzs6OunXr0rVrV1q1akWGDBkSFEBqooRJEl1kOGz9FjZ/DdERABiYmBNVn9ERHXiACyYTdK/qw3sNi+LqZG/jgEVEUp9k6/T94osvMmXKFO7cuUPZsmX5+uuvuXTpEqtXr6Zz586JniyNHDmSihUr4ubmRo4cOWjZsiWnTp164jbHjh2jTZs2+Pj4YDKZ+O6772KtN3HiRHx8fHB2dqZy5crs3r07UWMXeSb2jlB7KPTbArkrAGDCoJvdX2zMMJSa5kMYBszcfp4G4zaz+e+bNg5YRCR9eqqEKXfu3Lz33nscOXKE/fv388477+Dl5ZVkQfn5+dG/f3927tzJ2rVriYiIoEGDBjx48CDObYKDgylQoACjRo2KM7aFCxfyzjvvMHz4cPbv30/ZsmVp2LAhN27cSKpTEXk6OYrDa39Bw5HgYPkDJEf0TWY7jmac42Q8uM+VeyF0m76bwYsOcS843MYBi4ikL091S84wDJsOqnfz5k1y5MiBn58fNWvWjLe+j48PgwYNYtCgQTHKK1euTMWKFZkwYQIA0dHR5M2bl4EDBzJ06NB496tbcpIs7pyzPEl3brO16J45M0NDu7M6uhIA2TJapldpXFrTq4iIxCfZbsnZegTigADL9BJZsiS842t4eDj79u2jXr161jKz2Uy9evXYsWNHrNuEhYURGBgY4yWS5LLkh27LoNn34GT5YGeKvstkx+/42ek7snOXW/fDeH3efvrN2ceNwFAbBywikvY9VcLUvHlzvv/+++c60Pjx42nevPkzbxcdHc2gQYOoXr06pUqVSvDxb926RVRUFJ6enjHKPT098ff3j3WbkSNH4uHhYX3lzZs3wccXeSYmE5TvbhnwsmgTa3F90242ZRhCWzs/wGD1MX/qfevHor2XNEisiEgSeqqEacWKFRw8ePC5DnTw4EFWrlz5zNv179+fo0ePsmDBguc6fkIMGzaMgIAA6+vSpUvJHoOkc+65oON8aDsdMmQDwDX6Pl87/MQvzmPIY7pJYGgk7/12mO4z9nD5brCNAxYRSZue+hllf39/Nm/eHH/FJ2z/rAYMGMCKFSvYvHkzefLkSfCxAbJly4adnR3Xr1+PUX79+vU4O4k7OTnh5OT0XMcVeW4mE5RqA/lrw+ohcGQRAFU5xHrnIXwZ3oE5UfXZ/PdNGozbzJBGxehaxRuzWZP5iogklqdOmNasWcOaNWsSfKBn6ThuGAYDBw5k8eLFbNq0ifz58yf4uA85OjpSvnx51q9fT8uWLQHL7b7169czYMCA596/SJJzzQptpkKptrDibQi6ipMRymcOs2jtuJu3Q3txLjwnw5cdY/mhq4xuW4aC2TPaOmoRkTThqRKm7t27J3UcMfTv35/58+ezdOlS3NzcrK1THh4euLi4ANCtWzdy587NyJEjAUun7uPHj1v/feXKFQ4ePEjGjBkpVKgQAO+88w7du3enQoUKVKpUie+++44HDx7w6quvJuv5iTyXoo3Auyqs/QT2zQTA1zjBX87D+Dq8DVOjmrD3wl0aj9/CW3UL06dmARzs0ubo+yIiySVRJt9NbHG1RM2YMYMePXoAULt2bXx8fJg5cyYA58+fj7UlqlatWmzatMm6PGHCBMaOHYu/vz++vr58//33VK5c+ani0rACkuKc2wzLBsaYzPekuRBvhfTilJEPgJK53BnTtgwlc3nYKEgREduy6eS76ZESJkmRwh/Ahn8n88XycY4y2TMhojkTIlsSgT32ZhP9ahVkYF1N5isi6Y8SpmSmhElStEt7YGl/uPVoGqFzZm/eCunFYaMgYJnMd0zbMryQL7OtohQRSXbJNnCliKQCeSta5qSr+R6YLK1I+aMvsMRpOMMcfsGJcM7cuE+bSdv5bPlxgsMjbRywiEjqoYRJJC2xd4I6H0GfTeBVBgAz0fS1W876DB9SwXQSw4Dp287R6LstbD97y7bxioikEkqYRNKinGWg9wao8zHYOQKQJ/oKi5w+5zPH2WQglIt3gnnl5118sPgIQaERNg5YRCRlU8IkklbZOUDNwdBvK+SpCIAJg27m1WxyHUY181EA5u+6SINxm9l46oYtoxURSdGUMImkddmLQs810PArsLeMY5Yj6jrzHb9ijNM03AjmWkAor87Ywzu/HuRecLiNAxYRSXmUMImkB2Y7qNofXt8G3jWsxe1N69nkOpTa5gMA/LH/CvW+3czqo9dsFamISIr0XAnT8ePHefvtt6levTpFixbl/ffft67bvn0733//PXfu3HnuIEUkkWQtCN2XQ9NvwNEybUrWqFvMdBzLeOfJeHCfW/fD6Dd3P/3n7efW/TAbBywikjIkOGH69ttv8fX1Zfz48ezYsYMzZ85w61bMJ27efvttFi1a9NxBikgiMpuhYi94YwcUrGMtbsFm/FyH0sC8B4CVR65R/1s/lh68goZrE5H0LkEJ08qVKxk8eDB58+bljz/+4MaNG499oVarVo3s2bOzdOnSRAlURBJZpnzQ5Q9oPgGcLNOmZIq6wxTHcfzkPIEsBHI3OIK3Fhyk9+y9+AeE2jhgERHbSVDC9O233+Lq6sratWtp2bIl2bJli7Wer68vp06dinWdiKQAJhO80BX674QijazFDdmOn+sQXjbvAAzWnbhB/XF+/LrnklqbRCRdSlDCtG/fPqpUqUKBAgWeWC9btmz4+/snKDARSUbuuaDTAmj9M7hYpk1xiwpgguMPTHf5nmwEEBQayfu/H6bb9N1cuRdi44BFRJJXghKm8PBw3Nzc4q1348YN7O3tE3IIEUluJhOUaQ/9d0Px5tbiOsYu/FyH0MK8FTDYcvoWDb71Y+7OC0RHq7VJRNKHBCVM+fPn59ChQ0+sEx4ezuHDhylSpEiCAhMRG8mYAzrMgXYzIYPldrtrVCDjHX9kjss4cnCXB+FRfLTkKJ2n7uLi7WDbxisikgwSlDA1b96c8+fP8+2338ZZZ8yYMdy8eZPWrVsnODgRsaGSraD/LijVxlr0orEXvwxDaGvnBxjs+Oc2Db/bzMxt59TaJCJpmslIQA/Ou3fvUrZsWa5cuULbtm1p1aoVr7zyCo0bN6ZXr14sXryYefPmkT9/fg4cOPBUt+9Sg8DAQDw8PAgICMDd3d3W4YgknxMrYMXb8ODR9Ck7zC/wTvCrXCMrABV9MjOmbVnyZ3O1VZQiIrFKjN/vBCVMAH///Tdt27bl6NGjmEwmDMPAZDIBYBgGJUqUYMmSJRQqVChBgaVESpgkXQu+A6uHwuGF1qJQcwaGh3VmYVRtwISTvZnBDYrSs0Z+7Mwmm4UqIvL/bJowAURHR7N8+XL++usvzp8/T3R0NHny5KF+/fq0adMGOzu7hO46RVLCJAKcWg0rBkHQo+lTdpt9GRTck6tY+jyVy5eJsW3LUihHRhsFKSLyiM0TpvRGCZPIv0LuwZoP4OA8a1GYOQOfhnViflQdwISjvZm36xWh94v5sbfTtJUiYjtKmJKZEiaR/zi9Fpa9CUFXrUV77cry1oOeXCE7AGXzeDC2XVmKeKaNvowikvrYLGHavHnzU9VzdHQka9asFCpUyNq/KTVTwiQSi9AA+Osj2D/bWhRmzsDnYR2ZG1UXMOFoZ+ateoXpW7OAWptEJNnZLGEym83PlABlzJiRTp068eWXX5I1a9ZnPVyKoYRJ5AnOrINlb0HgZWvRfrsyvBn8GpcNS2tT6dwejG1XhmJe+vyISPKxWcLUo0cP7t27x7JlyzCbzfj6+pIvXz4ALl26xIEDBzAMg5dffpmQkBAOHjzIrVu3KFSoELt27SJz5swJCtbWlDCJxCM0AP76GPbPshaFm12srU0GZhzsTAysU5jXaxfEQa1NIpIMbJYw3bx5k8qVK1O0aFG+//57ChcuHGP9mTNnePPNNzl58iQ7d+7E1dWV3r17s3DhQoYMGcJXX32VoGBtTQmTyFM6s97St+n/WpsO2JVhYHBPLhs5ACiZy52v25WleE59lkQkadksYerVqxd//vknZ8+exdnZOdY6ISEhFCpUiEaNGjFt2jQCAwMpUKAAXl5eHD16NEHB2poSJpFnEBoIaz+GfTOtReFmFz4P78TcyDrW1qYBLxXmjZfU2iQiSScxfr8T9A21cuVKatasGWeyBODi4sKLL77IqlWrAHB3d6dcuXKcO3cuQYGKSCrj7A7NxkPXxeCeBwDH6BA+t5/OYtcx5DHdICLKYNy6v2kxYRvHrwbaOGARkbglKGEKCAggICAg3nqBgYEx6mXLli0hhxOR1KxgHXhjB5TvYS3yjTrMBpdhdLVfh4lojl8LpPmErXy37m/CI6NtF6uISBwSlDAVLlyYjRs3cuLEiTjrnDhxgo0bN1KkSBFr2bVr15Q0iaRH/9/a5JEXeNTa9IfrGPKYbhIZbfDdutO0mLiNY1fj/4NMRCQ5JShh6t+/P2FhYdSsWZORI0dy+vRpQkNDCQ0N5fTp04waNYpatWoRHh5O//79AUufpn379lG+fPlEPQERSUUK1oHXt8ML3a1F5f5tbepivx4wOHEtkBYTtjFurVqbRCTlSPBI32+//Tbjx4+PczwmwzAYNGgQ3377LQAnT55k/PjxtGnThnr16iU8YhtSp2+RRBTLk3T77X0ZeP9V6yjhxXO68027spTIpc+biCSczadG2bZtG5MmTWL79u1cu2aZiDNnzpxUr16dvn37UqNGjYTuOkVSwiSSyEID4a8PY44SbufKiLBX+CWyNmDC3mxiQJ1C9H+pkJ6kE5EEsXnClN4oYRJJImfW/dvadMVatNf+BQbef5VrWGYHKJHTMm6TWptE5FnZbFgBEZFEVaie5Um6cl2sRRUi97PJdSjt7f0Ag+PXAmkxcSs/rD9NRJT6NolI8lIL0zNQC5NIMji91tLaFHTVWrTLvgJv3u/BdbIAljnpvm5XlqJebraKUkRSEZu2MAUHB/PFF19QsWJFMmXKhJ2dXawve3v7hB5CRNKjwvUtrU2+na1FlSP34uc6jDZ2mwGDI1cCaPbDViZuPEOkWptEJBkkqIUpICCAF198kWPHjmFnZ4ejoyPBwcHkzJkTf39/Hu7S29sbIM2M7q0WJpFkdmo1LH8T7l+3Fm2zr8yg+925SSYAyuaxtDYV9lRrk4jEzmYtTKNGjeLo0aP06dOHwMBA2rZti8lk4sqVKzx48ICZM2fi5eVF5cqV+eeffxIUmIgIRRvBGzuhdHtrUfXIXWzOOJTmdtsBg0OXA2j6w1Ym+50lKlo9DEQkaSQoYVqyZAm5cuXi+++/x9nZOcZYTM7OznTr1o1169axePFivvnmm0QLVkTSoQxZoM3P0GEuZLDMFOASGcj3DhOYlXEiWQgkPDKaUX+epO3k7Zy9ed/GAYtIWpSghOnChQu88MILODg4WHZituwmIiLCWqdEiRLUqlWLmTNnPn+UIiLFm0H/XVCipbWoVuR2tmQcRkO7PQAcuHiPJuO3MHXLP2ptEpFElaCEydnZGWdnZ+vyw/uB/v7+MeplyZIlzfRfEpEUwDUbtJ8FbaeDS2ZLUeRdfnIYx5SMU3DnPmGR0Xyx8gQdp+zg/K0HNg5YRNKKBCVMefPm5dKlS9blYsWKAeDn52cti4yMZM+ePWTNmvU5QxQR+Y9SbeCNXVC0ibWoQeQmtrp9QG27gwDsOX+XxuO3MGv7eaLV2iQizylBCdOLL77I4cOHCQoKAqBZs2bY29vz5ptvMnnyZJYvX07btm05f/48tWrVStSARUQAcPOEjvOh5WRw8gDAPeIWMx3G8IPrDDISTEhEFMOXHaPLtF1cuhNs44BFJDVL0LACW7du5b333mPEiBE0bNgQgG+//ZbBgwdbO4AbhoGXlxd79+4lV65ciRu1jWhYAZEUKuAKLBsIZ9dbi+45ePL6g17siC4JgKujHR+9XIKOFfPGOWm4iKRNKW4uuV27drF48WLu3r1LkSJFePXVV8mSJUti7d7mlDCJpGCGAftnwZoPIfzRk3K/2zXmowftCMHS77JWkeyMalOanB4utopURJJZikuY0jolTCKpwN0LsLQ/nN9iLbrlmIc+Qb3YbxQBwM3Znk+bl6RVudxqbRJJB2w2cGXPnj2ZPn16vPVmzpxJz549E3IIEZGEyewN3ZZBo9Fgb2lFyhZ+md+dP+OzDL/iRDhBoZG88+sh+szZx82gMBsHLCKpQYISppkzZ7J169Z4623bto1Zs2Yl5BAiIglnNkOVftBvK+SpCIDJiKZb9BI2eXxKSdN5ANYev06DcX6sPHzNhsGKSGqQ4Ml3n0Z4eDh2dnZJeQgRkbhlKwQ910C9EWDnCEDOsHMsd/6EIS5LsSOKu8ER9J+/nzd/OcC94HDbxisiKVaSJUyGYbB//36yZ8+eVIcQEYmf2Q5qvA19NoFXaUuREcnrxkLWe3xBQdMVAJYdukr9cZvZcPL6E3YmIunVU3f6rlOnjvXfmzZtwsvLyzpg5X9FRkZy9uxZ/P396dq1a5qZHkWdvkVSuchw2DwWtnwDRhQAUWZHvo3uyI+hDTD+/RuyQ4W8fPRycdycHWwZrYgkkmR9Su7hfHEAJpOJ+DZzcHCgUaNGTJs2jWzZsiUouJRGCZNIGnF5HyzuC7dPW4tOOpWmV+BrXDZyAJA7kwtj25WhWsG08f0lkp4la8J04cIFwHKrrUCBArRt25axY8fGWtfR0ZFs2bJZJ+dNK5QwiaQhESGw/jPY+eOjIrsMfB7RhdnhtQDLcAM9qvkwpFExXBzVH1MktbLZOEyffvop5cqVo3nz5gk6aGqlhEkkDTq3BZa8AQEXrUX7nCrRL6A7N7FM8Fsgmytfty/LC/ky2ypKEXkOGrgymSlhEkmjQgNhzTA4MPdRkYMHQ0JfZWlEJQDMJnijdiHerFsYR/skfcBYRBKZzQauFBFJU5zdocVE6LQAXC19mJwjAhhv9x2zPKbgzn2iDZiw8QwtJ27jpH+gjQMWkeT2VAmTnZ1dgl/29vZJfQ4iIomjaGN4YyeUaGEtqhW2ie3uH/KS3WEAjl8LpPkP25jsd5aoaDXQi6QXT3VLzsfH57nmWzp37lyCt01JdEtOJJ0wDDjyG6x6F0IDrMXLHJswJLCtdSLfCt6Z+aZ9WbyzutoqUhF5CurDlMyUMImkMwFXYOkb8M8ma9Ed57y8FtiLA9GFAcjgaMeHTYvzSqV8mshXJIVSHyYRkaTkkRu6LIYmX1sn8s0Seok/nD7l04yLcSCS4PAoPlx8lB4z9nA9MNTGAYtIUlHCJCLyJGYzVOptmcg3dwXAMpFv98hFbMr8BYVMlwHw+/smDcZtZvmhq7aMVkSSyHMlTIcPH6Zv376UKFECDw8PPDw8KFGiBP369ePw4cOJFaOIiO09nMi3zkdgtjzMkjvkb9a4fMSbGdZiIpqAkAgG/nJAE/mKpEEJ7sM0fvx43nvvPaKiomKdJsXe3p6xY8fy1ltvPXeQKYX6MIkIAFcPWqZWuXnSWnTKpRyv3n2Vq1imUvF0d2Js27LULKIJyEVszWZ9mNauXcvbb7+No6Mjb7/9NgcOHODu3bvcu3ePgwcP8u677+Lk5MQ777zD+vXrn3n/I0eOpGLFiri5uZEjRw5atmzJqVOn4t1u0aJFFCtWDGdnZ0qXLs2qVatirO/RowcmkynGq1GjRs8cn4ikc7l8oc8mqNLfWlQ05AB+GT+gk/N2wOB6YBjdpu/m4yVHCQmPslWkIpJIEtTC1LhxY9avX8+mTZuoVq1arHV27NhBzZo1qV+//mOJS3waNWpEx44dqVixIpGRkXzwwQccPXqU48eP4+oa++O727dvp2bNmowcOZKXX36Z+fPnM3r0aPbv30+pUqUAS8J0/fp1ZsyYYd3OycmJzJmfbroDtTCJyGP+8bNMrRJ42Vq0y6Umfe925h5ugGVqlW/al6WcplYRsQmbDSuQNWtWXnjhBdauXfvEevXr12f//v3cvn07QcE9dPPmTXLkyIGfnx81a9aMtU6HDh148OABK1assJZVqVIFX19fJk+eDFgSpnv37rFkyZIExaGESURiFXIP/nwfDi+0FgU7ZefN4NdYF1EGADuzif61CzKwbmEc7PS8jUhystktueDgYLJnj/++fPbs2QkODk7IIWIICLAMHJclS5Y46+zYsYN69erFKGvYsCE7duyIUbZp0yZy5MhB0aJFef3115+YzIWFhREYGBjjJSLyGJdM0HoKtJsJLpZWpAxhN5lqN4ofM83DmTCiog2+33CG1j9u58yNIJuGKyLPLkEJU968edmxYweRkZFx1omMjGTHjh3kzZs3wcEBREdHM2jQIKpXr269tRYbf39/PD09Y5R5enri7+9vXW7UqBGzZ89m/fr1jB49Gj8/Pxo3bkxUVOz9C0aOHGl9+s/Dw+O5z0VE0riSreD1HVCwrrWoSehKtmcejq/dPwAcuRJA0++3Mmv7eaI1tYpIqpGghKlFixZcuHCBnj17cu/evcfWBwYG0rt3by5evEjLli2fK8D+/ftz9OhRFixY8Fz7AejYsSPNmzendOnStGzZkhUrVrBnzx42bdoUa/1hw4YREBBgfV26dOm5YxCRNM49J3T5PeZglyEXWew4gk/cV2BHFGGR0QxfdozuM3bjH6DBLkVSgwTNjDts2DD++OMP5s2bx9KlS2nUqBE+Pj4AXLhwgdWrVxMYGEiBAgUYNmxYgoMbMGAAK1asYPPmzeTJk+eJdb28vLh+/XqMsuvXr+Pl5RXnNgUKFCBbtmycOXOGunXrPrbeyckJJyenhAUvIumXyWQZ7LJAbfijN1w9gMmIpGf4fOpnOUSXu69xwfBiy+lbNPxuM1+2KsXLZXLZOmoReYIEj8N09epV+vbty8qVK2Nd37RpU3766Sdy5Xr2LwHDMBg4cCCLFy9m06ZNFC5cON5tOnToQHBwMMuXL7eWVatWjTJlylg7ff/X5cuXyZcvH0uWLKF58+bxHkOdvkXkmUVFgN8Y2PI1GNGWIvsMjDK68/ODGoBl/rmWvrn4tEUpPFwcbBisSNqUIibfPXfuHFu3buXqVct0ALly5aJGjRrkz58/wft84403mD9/PkuXLqVo0aLWcg8PD1xcLE3c3bp1I3fu3IwcORKwDCtQq1YtRo0aRdOmTVmwYAFfffWVdViB+/fv8+mnn9KmTRu8vLw4e/Ys77//PkFBQRw5cuSpWpKUMIlIgl3aY2ltunvOWnTYtRo9bnfjDpbvk1weznzdvizVCmazVZQiaVKKSJiSQlwzfs+YMYMePXoAULt2bXx8fJg5c6Z1/aJFi/joo484f/48hQsXZsyYMTRp0gSAkJAQWrZsyYEDB7h37x65cuWiQYMGfP755491Fo+LEiYReS5h92HNB7B/lrUo1Ckr74T1ZlWoZfgBkwl61cjPuw2K4uxgZ6tIRdIUmyVMgwcPpmvXrpQtWzZBB02tlDCJSKI4uQqWDYTgW9aivzK8zJt32hCKpbW7mJcb33X0pZiXvmtEnpfNEiaz2YzJZKJ48eJ07tyZV155BW9v7wQFkJooYRKRRHP/BiwdAKfXWIvuZfDh1YA+HIjyAcDRzsz7jYrSs3p+zObYW95FJH42S5h++OEH5s2bx+7duy07MZmoXr06nTt3pn379k891Uhqo4RJRBKVYcC+GbD6A4gMsRSZ7Znp2InP7zUk+t+RX6oVzMo37cuS08PFltGKpFo278P0zz//MHfuXObNm8fp06cxmUw4ODjQqFEjOnfuTPPmzdPUY/lKmEQkSdw6bR1+4KGLGcvyyu2eXDYssyq4O9vzZavSNCur4QdEnpXNE6b/t2/fPubOncvChQvx9/fHZDLh5uZG69atmT59emIcwuaUMIlIkomKgE2jYOu31uEHIh0y8rnxGrPuV7ZWa10uNyNalMTdWcMPiDytFJUwPRQdHc2GDRuYPn06CxYswGQyxTn1SGqjhElEktyFHbC4D9y7aC3a61aXnjc7EogrALkzuTCugy+V8sc9v6aIPGKzyXefZPPmzfz666+sWbMm/soiIhKTd1XotxXKdLAWVQhaz87MH/OS098AXLkXQocpOxi9+iThkdG2ilQkXUmUhOngwYO899575M2bl7p16zJ16lQiIyPp1q0bq1evToxDiIikH84e0HoKtJkGTh4AZAjxZ7rpU8ZlXYoDkRgGTNp0ljaTtnP25n0bByyS9iX4lty5c+eYP38+8+fP5+TJkxiGgYODAw0bNqRz5860aNECZ2fnxI7XpnRLTkSS3b1LsLgfXNhqLbrhVoJXbvfiTLRlrkxnBzMfNS1B58r54hz4VyQ9s1kfpqpVq7J7924eblqtWjU6d+5Mhw4dyJIl7d5TV8IkIjYRHQXbv4cNX0B0pKXI3oVxdj35IaAaD+ejq1c8B6PblCFrxrTzdLJIYrDpwJXFihWjc+fOdO7cGR8fnwQdPLVRwiQiNnX1APzeG26fthYd86hJ5+uduYcbANkyOvF1uzLULprDVlGKpDg2S5gOHDhAuXLlEnTA1EwJk4jYXPgDWPOhZcDLf4W65GBQWD9WBxezlvWo5sPQxsU0H50IKXRYgbRMCZOIpBgnVljmowu5A4CBiT/d2/HWjZeJwB6Aop5ujO+k+ehEUuSwAiIikgyKvwyvb4cCtQEwYdAk8Fd2Zh9JUXt/AE5dD6L5hG1M33oO/W0s8nyUMImIpFbuOaHLYmjwBZgtI39nDTrBn84fMijzdsAgPDKaz1Ycp8eMPdwICrVtvCKpmBImEZHUzGyGagOh93rIWthSFBnCoJAJ/JlzKh5Yxmjy+/smjb/bwvoT120ZrUiqpYRJRCQtyFkW+vpB+R7WouJ3N7Iryyc0dD0DwO0H4bw2ay/Dlx4lNCJtTFklklyUMImIpBWOrtBsPHSYCy6ZAXAO9mdy1HB+8FqJPZYxnGbtuECLCds45R9ky2hFUhUlTCIiaU3xZtBvG/i8CFg6hDe7N48dnl9T0P4mYOkQ3mzCVmZtP68O4SJPQQmTiEha5JEbui2FusPBbBlmIHvAYf5y+ZB+WfYBEB4ZzfBlx3ht1l5u3Q+zZbQiKd5TjcM0e/bs5zpIt27dnmv7lELjMIlIqnR5H/zeE+6etxYdzNqEzlfa8AAXALK7OfFt+7K8WDi7jYIUSTrJNnCl2WxO0ISOhmFgMpmIikobnQuVMIlIqhUaCKsGw+GF1qLgjN70CXmDrQ/yWsv61CzA4AZFcbTXDQhJOxLj99v+aSp98sknmgFbRCQ1c3aH1lOgYF1Y+Q6E3yfD/QvMMX/EAq9X+cC/JgZmpmz+hx1nbzO+oy8Fsme0ddQiKYamRnkGamESkTThzj/w22twdb+16FLWarTz74p/lAcAGRzt+LR5SdqWz6M/mCXV09QoIiLy7LIUgJ5roPogwJIM5b29nS3un9A+898ABIdH8d5vh3lrwUECQyNsF6tICqGESUQkPbJ3hPqfQtc/IKMnAA4hNxkTMoKZeVbg8O+YTcsOXaXp91vYf/GuLaMVsbnnuiV38eJFli9fzunTpwkKCop1LA+TycS0adOeK8iUQrfkRCRNun8TlvSDM+usRXczl+GVu705EZoVADuziXcbFKFfzYKYzbpFJ6lLsj0lF5vPPvuMzz//nOjoaGvZw109vN+tp+RERFKJ6GjY+SOsGwHRlltw0Y5ujHN+gx9ulLVWq14oK+Pa+5LD3dlGgYo8O5v1YVq4cCEjRowgb968TJkyhfr16wOwZs0aJk2aRK1atTAMg3feeYcNGzYkKDAREUlGZjNUGwC91kLm/Jai8CDeDRzNsnwLcTFZBrbcduY2jcdvYeOpG7aMViTZJShh+vHHH3F0dGTjxo289tpr5MyZE4D69evTt29fNmzYwDfffMP48eOxs7NL1IBFRCQJ5SoHfTdD6XbWojI3lrIvx5dUzXgdsEzi++qMPXyx4jjhkdFx7UkkTUlQwnT48GGqVauGt7c3EPMW3ENvv/02RYsW5YsvvkiEMEVEJNk4u0Prn6HFj+CQAYAMAWeYzzA+zb0bsHzXT916jraTt3P+1gMbBiuSPBKUMIWFheHl5WVddna23Mu+d+9ejHply5Zlz549CY9ORERsw2SCcp0trU2epS1FkaF0v/0dG71nksUuFIDDlwN4+YetLD14xZbRiiS5BCVMOXPm5MaNR/evc+fODcCxY8di1Lt8+XKa6fAtIpIuZSsMvdZBxd7WovzX17Izy3AaZ7YkSffDInlrwUHeW3SI4PBIW0UqkqQSlDCVLl2aU6dOWZdr166NYRgMHz6cBw8sTbO//vorW7ZsoWTJkokTqYiI2IaDMzT9GtrPAWfLSOCOQZf4MWwY47238vAW3aJ9l2k+YRsn/QNtGKxI0khQwtSsWTOuXLlifQKuevXqvPTSS2zcuJHMmTOTLVs2OnXqhMlk4uOPP07UgEVExEZKNId+WyFPJQBM0ZG0uP4jO7ynkNMxGIAzN+7TYsI25u26EOvYfCKpVYISpi5dunDixAl8fX2tZYsXL6ZPnz5kyZKFoKAgSpQowZw5c2jUqFFixSoiIraWKR+8ugpqvG0tynndjy1uH9Eu20UAwiKj+XDxUQbMP6BpVSTN0OS7z0ADV4qI/J/T62BxXwi+BYBhMrPeqxe9z9XE+Pfv8bxZXJjQ6QXK5s1kw0AlvdPkuyIiYjuF61lu0fm8CIDJiKbetSns8Z6Mj7OlP+ulOyG0nbydqVv+0S06SdWUMImISMK554RuS6HWUMAyJl+261tZ7/oRXbwuARARZfDFyhP0nr2Xuw/CbRisSMIlOGE6fvw4PXr0oECBAri4uGBnZxfry97ePjHjFRGRlMZsBy8NsyROrjkAsHtwnc8DhjG70CbMWEYDX3fiBk2+38Le83dsGa1IgiSoD9OOHTuoV68eISEhAGTJkgU3N7c46587dy7hEaYg6sMkIhKPoOvwR28452ctuuNVnXY3enI22AUAO7OJwQ2K0rdmAcxmk60ilXQkMX6/E5Qw1a5dm82bNzNo0CA++ugjsmTJkqCDpzZKmEREnkJ0FGz5BjaNBMPSuhTl6slnTu8y62oea7XaRbPzbXtfsrg62ipSSSdsljBlzJiRIkWKsH///gQdNLVSwiQi8gzObYHfX4P7lkl7DZOZbXn70O10DaINS48QL3dnvu9Ujkr508cf3mIbNntKztHRkWLFiiXogCIikk7kf9HyFF2B2oDlKboaFyez1+cnCmawzEXnHxhKp5938uOmM0RH6yk6SbkSlDDVqFGD48ePJ3YsIiKS1mTMAV3+gJc+BJPlJyfLtS38leFDuuW+BkBUtMGY1ad4bdYe7ugpOkmhEpQwffXVV5w5c4aJEycmdjwiIpLWmO2g1vvQdcmjp+juX+PTO+8zq9hOTCZLy9LGUzdpqqfoJIVKUB+m2bNns2fPHn788Udq1KhB/fr1yZMnD2Zz7PlXt27dnjvQlEB9mEREnlPQdUu/pvNbrEW3c9ehzbVunA+2dP62M5sY0qgovV8sgMmkp+jk+dms07fZbMZkMsUYtTW2N7VhGJhMJqKiohIUXEqjhElEJBFERYLfKNg89lGRe14+cHyfhZezWsvqFffkm3Zl8cjgYIsoJQ2xWcI0YsSIZ8r6hw8f/qyHSJGUMImIJKLT6yxjNoVYbsEZdo6syzeI3ifK8nDU8DyZXZj4iuaik+djs4QpvVLCJCKSyAIuw6IecHmPtei698u0utieqyGWmSIc7Ex81LQE3ap66xadJIgm3xURkdTNIw/0WAVV3rAWeV5YgV/mz2mWKxCwzEU3fNkxBv5ygPthkbaKVNI5JUwiImJb9o7QaCS0nw1Olr/+He6c5vugdxlX/LS12orD12g+YSun/INsFamkY091S65nz56YTCa++uorPD096dmz59MfwGRi2rRpzxVkSqFbciIiSez2Wfi1G1w/ai26ULAzrc404U6Y5Xacs4OZr1qVpvULeeLai0gMydaH6eFTcSdOnKBIkSJxDh8Q6wH0lJyIiDyL8GBYNRgOzrMWhXm+QO/QgWy+7mQt61QpH8OblcDZwc4WUUoqkhi/3/ZPU2njxo0A5MuXL8ayiIhIonPMAC0mQt7KsOo9iArD6fp+Zrm8y7RiH/HFSS8Aftl9kSNX7jGpc3nyZslg46AlrdNTcs9ALUwiIsns6kHLLbp7F/4tMHG06ADaHqtK6L/9v92d7fmuoy91innaKkpJ4fSUnIiIpG25fKGvHxRu+G+BQalTP7C3wFRKZ7F09wgMjaTnzL2MXXOSKE3gK0lECZOIiKRsLpmh0wKo8xEPB7TMeHEDSx0+olehQGu1iRvP0m36Lm7fD7NRoJKWJThhOnbsGK+++ioFChTAxcUFOzu7WF/29k/VTUpERCRuZjPUfA+6/gEZLNOnmAMu8KH/IGaWO42d2ZJIbTtzm5d/2MqBi3dtGa2kQQnqw+Tn50fjxo0JDQ3FZDKRJUsWMmbMGGf9c+fOPVeQKYX6MImIpAD3Lln6NV3dby26XuQVWv3TnKv3owHL6OCfvFyCLlU0OrjYcGqUypUrs2fPHj766CMGDx6cbpIHJUwiIilEZBj8OQT2zbAWRXi9wICoQay59OjORqtyufmqVWlcHDX0QHpms4TJxcWFcuXKsX379gQdNLVSwiQiksIcmAsr3oEoS78lI0M25ub5lI8PZ7ZWKeblxk9dy+Od1dVWUYqN2ewpuaxZs+Lj45OgAz6NkSNHUrFiRdzc3MiRIwctW7bk1KlT8W63aNEiihUrhrOzM6VLl2bVqlUx1huGwSeffELOnDlxcXGhXr16nD59Oo69iYhIileuC7z2F2SyjBNoCr5F19NvsqriYVwdLT9xJ/2DePmHraw/cd2WkUoql6CEqWnTpuzcuTPJRvD28/Ojf//+7Ny5k7Vr1xIREUGDBg148OBBnNts376dTp068dprr3HgwAFatmxJy5YtOXr00fD6Y8aM4fvvv2fy5Mns2rULV1dXGjZsSGhoaJKch4iIJINcvtDHDwrWsSwbUZQ4MoqdReZTIpvlVlxQaCSvzdrLt3+d0tADkiAJuiV38+ZNqlSpwosvvsj48ePx8PBIithiHC9Hjhz4+flRs2bNWOt06NCBBw8esGLFCmtZlSpV8PX1ZfLkyRiGQa5cuXj33XcZPHgwAAEBAXh6ejJz5kw6duwYbxy6JScikoJFR8HGL2HLN9aiqOzF+dTlA2b//agPU+2i2RnfoRweGRxsEaXYQLJNjfJf2bNnZ/fu3dSqVQsfHx8qVKhA7ty5Y51jLjEm3w0ICAAgS5YscdbZsWMH77zzToyyhg0bsmTJEsDypJ6/vz/16tWzrvfw8KBy5crs2LEj1oQpLCyMsLBH43kEBgY+VkdERFIIsx3U/QRyl4c/+kJ4EHY3T/Cp8wAqV/qCgXuyEm3AplM3aTZhKz91LU/xnPrjV55OghKmwMBA2rVrx4kTJzAMg/Xr18dZ93kTpujoaAYNGkT16tUpVapUnPX8/f3x9Iw5LL6npyf+/v7W9Q/L4qrzXyNHjuTTTz9NcOwiImIDxZpCn42w4BW49Tem0ACaHn6TUpXeptXhytwJieLinWBa/biN0W3K0MI3t60jllQgQQnT4MGD2bRpE6VKlaJ3794UKFDgieMwPY/+/ftz9OhRtm7dmiT7f5Jhw4bFaLUKDAwkb968yR6HiIg8o2yFodd6WPI6nFwBGHgf+pZtBRvR/c5r7L4WQWhENG8tOMjhywEMa1wMeztNfiFxS1DCtHTpUvLmzcuOHTtwdU26xzQHDBjAihUr2Lx5M3ny5HliXS8vL65fj/kExPXr1/Hy8rKuf1iWM2fOGHV8fX1j3aeTkxNOTk7PcQYiImIzzu7Qfg5s/RY2fAEYuJxdzYJs/zC21HAmHbUMaDlt6zmOXQ1g4isvkDWjvvMldglKp0NCQqhSpUqSJUuGYTBgwAAWL17Mhg0byJ8/f7zbVK1a9bFbg2vXrqVq1aoA5M+fHy8vrxh1AgMD2bVrl7WOiIikMWYz1BwMnReBs+UBJfOtv3n/0utMr34XBztL0rTznzs0+2ErRy4H2DJaScESlDD5+vrG2e8nMfTv35+5c+cyf/583Nzc8Pf3x9/fn5CQEGudbt26MWzYMOvyW2+9xerVq/nmm284efIkI0aMYO/evQwYMACw9KUaNGgQX3zxBcuWLePIkSN069aNXLly0bJlyyQ7FxERSQEK14feGyF7cQBMYYHU2TeATZX3kT2jIwBXA0JpM3k7v+27bMtIJYVKUML0ySefsH37dlavXp3Y8QAwadIkAgICqF27Njlz5rS+Fi5caK1z8eJFrl27Zl2uVq0a8+fPZ8qUKZQtW5bffvuNJUuWxOgo/v777zNw4ED69OlDxYoVuX//PqtXr8bZ2TlJzkNERFKQrAWh11oo3uzfAoPc+79mS/5ZVMtjuRUXHhnN4EWHGLHsGBFR0baLVVKcBI3DtHnzZv744w8mTpxI586dqV+/fpzDCgBxjp2U2mgcJhGRNMAwYMvXsOFLwPITGJ2jJN9kGc7Eg5HWalULZGXCK+XUrykNsNlccmazGZPJxMNN45sJOqlGBE9uSphERNKQv9fA770g7N8x9lyysKH0GPpucyUiyvL7ljuTC1O6ladkrqQdoFmSls0Sph49esSbJP2/GTNmxF8pFVDCJCKSxtw6Db90hNtnLMsmOy5W+pA2+8pw8344AM4OZsa0LUvzsrlsGKg8D5slTOmVEiYRkTQo5B780RtO/2UtCi7Zia7XO7Lv8qM5TF+vXZDBDYpiZ376BgNJGRLj9ztBnb5feOEF2rVrl6ADioiIpCgumaDTAqjxaKDiDMd+YZHzl/Qs42Itm7TpLL1m7SEwNMIGQYqtJShhOnXqFA4OmrRQRETSCLMd1BsObaaBveXJafPl3XzsP4AfapusrUobT92k5cRtnL1535bRig0kKGEqXLgwt2/fTuxYREREbKt0W+i5Gtwt88uZAq/QbO+rrKpznUwZLA0F/9x8QMsJ29h48oYtI5VklqCE6bXXXsPPz4+TJ08mdjwiIiK2laucZZDLvJUty5GhFN06CL/y2yjuaZnhIigskp6z9jBl81nUFTh9SFDCNHDgQHr06EGtWrUYN24cZ86cITw8PLFjExERsQ03T+i+HHy7WIs89nzHcq+faV7cMsSAYcBXq07y7qJDhEakjeFzJG4JekrOzs4OsMz5Ft/wAiaTicjIyCfWSS30lJyISDpjGLDzR/jrIzAsI38bnqWYlnckX2wNslYrly8TP3UpTw53zRyREtlsWAEfH59nGofp3Llzz3qIFEkJk4hIOnV6LfzW89Egl67Z2Vbhe17bYCI0wpJIebk7M7V7BUrl1iCXKY3GYUpmSphERNKxm6dgfge4+28jgJ0Tl2qOocP2vFwNCAUsg1yOa+9L49I5bRio/JfNxmESERFJd7IXhd4bwOdFy3JUGHk3vsVa3y1UyGdpVQqNiOb1efsZv+60OoOnMYmSMIWFhXHt2jXu3LmTGLsTERFJmTJkga6L4YXu1iLXXeNYmGUyHX2zWMvGrfubgb8cUGfwNOS5EqYpU6ZQrlw5XF1dyZMnD4MHD7au++OPP2jdujVnzpx57iBFRERSDDsHaDYeGo0Ck+Vn1O7kckYGDOWzOll52MV3xeFrdPhpBzcCQ20YrCSWBCVMUVFRtGrVitdff50TJ05QvHjxx5oey5Yty5IlS1i4cGGiBCoiIpJimExQ5XV45VdwdLMUXTtItyOvsuDlDLg6Wp4mP3Q5gBYTt3HsaoAto5VEkKCEacKECSxdupTGjRtz4cIFjhw58lidggULUqhQIf7888/nDlJERCRFKlwfeq2FTPksy0FXqbypM6sb3yd3Jss8dNcCQmk3eQd/HfO3YaDyvBKUMM2cORNPT08WLlyIp6dnnPVKlCjBhQsXEhyciIhIipejOPTaAHkqWZYjHpB3TS/WVD5MubyWzuDB4VH0nbuPn/w0MnhqleDJdytXroyrq+sT67m6unLz5s0EBSYiIpJqZMxuGRm8VNt/Cwwy+g1nUZ5FtCyTw1JiwMg/TzLsjyNEREXbLlZJkAQlTA4ODoSGxt+J7eLFi7i5uSXkECIiIqmLgzO0mQq1hlqL7A/MZFzUSIbUzmUtW7DnEj1m7CYgOMIWUUoCJShhKlmyJPv27SMoKCjOOjdu3ODgwYP4+vomNDYREZHUxWSCl4ZBqylg52gpOrue1//pz5TmOXC0t/zsbjtzm9aTtnHh9gNbRivPIEEJU9euXbl9+zb9+vWLddLdqKgo+vfvT3BwMN27d49lDyIiImlY2Q7QdQm4ZLYs3zhGg+2dWdrKhSyulkTq7M0HtPpxO/suaAzD1OCpEqYCBQowZMgQ63KfPn2oXbs2v/zyC0WLFqVfv34AHDp0iLfeeosiRYrw+++/U79+fTp37pw0kYuIiKRkPtXhtXWQpYBl+f51iq/uxJrGDyiUIyMAdx6E0+nnXSw/dNWGgcrTeKqE6fz58zE6b9vZ2bFq1Spef/11rl69ypQpUwA4cOAAP/zwAxcvXqR3794sWbLkmSbpFRERSVOyFbIkTfmqWpYjgsm+8lWWVz5B9UJZAQiPjGbgLweYuPGMnqBLwZ5q8l2z2UyPHj2YPn36Y+tu3rzJpk2bOH/+PNHR0eTJk4eXXnqJXLlyxbKn1E2T74qISIJEhMLSN+Do79aiqCoD+CCoLQv3PWpd6lAhL1+0KoWDnaZ6TUyJ8ftt/7xBZM+enXbt2j3vbkRERNIuB2doPdUywOXWcQDY7ZzAqOIXyV9vCKPWWcYsXLj3ElfuhfBjlxdwd3awZcTyH0phRUREkoPZDPVGwMvfgckydYrpxDL6nX+bSa28cfy3VWnrmVu0n7yDawEhtotVHvPUt+R8fX1p2bJlgg7yySefJGi7lEa35EREJFGcXguLekD4fcty1kIcrj2NbotvcO/f8Zm83J2Z3qMiJXLp9+Z5Jcbv91MnTAnpvG0YBiaTiaioqAQFl9IoYRIRkURz7RDMawf3r1uWXbNzucksXlkZzsU7wQBkdLLnx84vULNIdhsGmvola8JUqFAhqlevnqCDzJgxI0HbpTRKmEREJFHduwhz28KtU5ZlhwwEvPwz3bdm5uClewDYm02MbF2adhXy2i7OVC5ZE6a4npJLT5QwiYhIogu+Aws6w8XtlmWTHeGNv6X/yVKsPX7dWu3d+kUYUKeQhutJgMT4/VanbxEREVvKkAW6LoYSLS3LRhSOq97ip3wb6FHV21rtm7V/88Hio0Rq4l6bUMIkIiJiaw7O0HYGVHnDWmTe9CXDzdP4oFFha9kvuy/Sb+4+gsMjbRFluqaESUREJCUwm6HhV1D/c2uRad90+vh/yg/tiuFgZ7kVt+7EDV75eRd3Hjw+l6skHSVMIiIiKYXJBNXfhNY/g/nfsaVPrqDZwTeY+0pR3JwsZQcv3aPtpO1c+vdpOkl6T9XpWyzU6VtERJLN2Q2wsOujsZpylOB0g5l0+fUS1wPDAMju5sSsVytprKZ4qNO3iIhIWlWwDvRYCa7/jsF04ziFl7dlacccFMzuCsDNoDA6/LSD7Wdv2TDQ9EEJk4iISEqVyxd6roHMPpblgIt4LWrB4hbOlMuXCYCgsEh6TN/DysPXbBVluqCESUREJCXLWhB6/gWepS3LIXdwX9iaBXVDqFssBwDhUdEM+GU/c3ZesGGgaZsSJhERkZTOzRNeXQneNSzLEQ9wWtiRKeUv075CHgAMAz5ecpTx606j7smJTwmTiIhIauDsAV1+h2IvW5ajI7D7/VVG+xygX62C1mrj1v3N8GXHiI5W0pSYlDCJiIikFg7O0G4W+Hb5t8DAtOIthmZcxYeNi1mrzd5xgTcXHCA8UqOCJxYlTCIiIqmJnT20mABVBzwqW/8pvcNm8E3bMtiZLQNcrjh8jd6z9xISHmWjQNMWJUwiIiKpjckEDb6AusMflW3/gTbXvmZKZ1+c7C0/735/36TrtF0EhETYKNC0QwmTiIhIamQywYvvwMvfAZZWJfbNpO6Jj5jTo5x1VPC9F+7SccpObgaF2SzUtEAJk4iISGpW4VVoM/XRVCpHf6fSrjdZ0NOXrK6OAJy4Fki7yZpK5XkoYRIREUntSreFjvPB3tmyfHoNJTf05LeepcnlYSk7fzuYdpN3cObGfRsGmnopYRIREUkLijSEzr+BY0bL8oWt5F/1Cr+/WoIC/06l4h8YSoefdnD8aqANA02dlDCJiIikFflfhO7LwCWzZfnKPnIubsdvXQpRIqdl0tnbD8LpOGUH+y/etWGgqY8SJhERkbQkd/l/J+21TJvC9SNkWdSKBZ28eeHf+ecCQyPpMnUX289o0t6npYRJREQkrfEsCa/+Ce65Lcu3/sb9l2bMbZOT6oWyAhAcHkWPmXvYcPK6DQNNPZQwiYiIpEXZClmSpkzeluW758kwrxnTm2WlXnFPAMIjo+k7Zx+rj16zYaCpgxImERGRtCqzN/RcDdmKWJYDL+M0txmTGrnRrGwuACKiDPrPP8DSg1dsGGjKp4RJREQkLXPPZenTlKOEZTnoGg6zX+a7l5xoWz4PAFHRBoMWHuTXPZdsGGjKpoRJREQkrcuYA7qvAK/SluUHN7Cb/TJjqpvoXDkfAIYB7/9+mDk7ztsuzhRMCZOIiEh64JoVui2DXOUsy8G3Mc9pzheVI+lZPb+12sdLjzF96zkbBZlyKWESERFJLzJkgW5LIU9Fy3LIXUyzW/Jx+TDeqF3QWu2zFceZuuUfGwWZMilhEhERSU+cPaDrYshX1bIceg/T7Ba8VzqEt+oWtlb7YuUJJvudtVGQKY8SJhERkfTGyQ06L/q/pCkA05wWvF3iPu/UL2KtNurPk0zceMZGQaYsSphERETSIyc3y9xz+apZlkMDYHZL3iwWyHsNi1qrjV1ziu/Xn7ZRkCmHEiYREZH0yimjpaXJu7plOSwAZreif5FAhjUuZq327dq/+SGdJ01KmERERNIzp4zwyq8xk6Y5rehb5D4fNS1urfbN2r+ZsCH9Jk1KmERERNK7hy1N1ttz92B2C3oVfsCHTR4lTV//9Xe67dOUIhOmzZs306xZM3LlyoXJZGLJkiXxbjNx4kSKFy+Oi4sLRYsWZfbs2THWz5w5E5PJFOPl7OycRGcgIiKSyji6QudfIW8Vy3LIXZjdgt5FQ/mgyaPbc2PXnEqXSVOKTJgePHhA2bJlmThx4lPVnzRpEsOGDWPEiBEcO3aMTz/9lP79+7N8+fIY9dzd3bl27Zr1deHChaQIX0REJHV6+PTcw3Gagm/DrGb0KR4Zo0/T2DWnmLI5fQ05YG/rAGLTuHFjGjdu/NT158yZQ9++fenQoQMABQoUYM+ePYwePZpmzZpZ65lMJry8vBI9XhERkTTD2R26/A6zW8LV/RB8C2Y1p++rqzAoxqg/TwLw1aqTONiZefX/RglPy1JkC9OzCgsLe+z2mouLC7t37yYiIsJadv/+fby9vcmbNy8tWrTg2LFj8e43MDAwxktERCTNezi4Zc6yluX7/jC7Bf3KOjC4waNxmj5dfpy5O9PH3Zo0kTA1bNiQqVOnsm/fPgzDYO/evUydOpWIiAhu3boFQNGiRZk+fTpLly5l7ty5REdHU61aNS5fvhznfkeOHImHh4f1lTdv3uQ6JREREdtyyQRdl0COEpblgEswqxkDKrjyZp1C1mofLTnKr3su2STE5GQyDMOwdRBPYjKZWLx4MS1btoyzTkhICP3792fOnDkYhoGnpyddunRhzJgx+Pv74+np+dg2ERERFC9enE6dOvH555/Hut+wsDDCwsKsy4GBgeTNm5eAgADc3d2f+9xERERSvPs3YEYTuP3vkALZimD0WMmoLbf5yc8y35zJBN+2L0urcnlsGGjcAgMD8fDweK7f7zTRwuTi4sL06dMJDg7m/PnzXLx4ER8fH9zc3MiePXus2zg4OFCuXDnOnIm7p7+TkxPu7u4xXiIiIulKxhzQfRlk9rEs3/ob0+yWDK3lSc9/+y8ZBgxedJjVR/1tF2cSSxMJ00MODg7kyZMHOzs7FixYwMsvv4zZHPspRkVFceTIEXLmzJnMUYqIiKQy7rmg2zJwz21ZvnEM07x2fFw/D12q5AMgKtpg4C/72XTqhg0DTTopMmG6f/8+Bw8e5ODBgwCcO3eOgwcPcvHiRQCGDRtGt27drPX//vtv5s6dy+nTp9m9ezcdO3bk6NGjfPXVV9Y6n332GX/99Rf//PMP+/fvp0uXLly4cIFevXol67mJiIikSpm9oftycM1hWb6yF9PCznzWpBCtX7AkUhFRBn3n7GPXP7dtGGjSSJEJ0969eylXrhzlypUD4J133qFcuXJ88sknAFy7ds2aPIGlteibb76hbNmy1K9fn9DQULZv346Pj4+1zt27d+nduzfFixenSZMmBAYGsn37dkqUKJGs5yYiIpJqZS1oeXrO2cOyfG4z5t97MqZlcZqUtgzbExYZzWuz9nLw0j3bxZkEUnyn75QkMTqNiYiIpHqX9sDsFhDxwLJcuj3hzSfRd+5+Np66CYCHiwML+1ahmJftfy/V6VtERESSX96K0OkXsHOyLB/5Fcc17zOp8wtULZAVgICQCLpO283F28E2DDTxKGESERGRZ1egFrSbCSY7y/LeaThvG8PU7hXwzZsJgJtBYXSetpPrgaE2CzOxKGESERGRhCnWBFpOerTsNxrXg9OY0aMiRTwzAnDpTgjdpu3mXnC4jYJMHEqYREREJOHKdoCGIx8t//k+mf9ZxpzXKpMnswsAp64H8erMPTwIi7RRkM9PCZOIiIg8n6pvwIvvPlpe3BfP61uZ16sy2d0s/ZwOXLxHv7n7CI+MtlGQz0cJk4iIiDy/Oh9D+R6Wf0dHwq9d8Q45weyelXB3tgdgy+lbDF50iOjo1PeAvhImEREReX4mEzT9Foo3tyxHBMP89hR3uMH0HhVxsrekHMsOXeXzlcdJbaMaKWESERGRxGG2gzZTwedFy3LwbZjbmgpZI5j4ygvYmU0AzNh2nsn/TtybWihhEhERkcRj7wQd50GOkpblexdgfjvqFczAyFalrdVGrz7Jr3sv2SjIZ6eESURERBKXswd0+Q088lqWrx2ChV1pX86T9xoWtVYb9scRNpy8bqMgn40SJhEREUl87rmgy+/gnMmy/M9GWDaQN2oVoEc1HwCiog36zzvA4cv3bBXlU1PCJCIiIkkje1F4ZSHYO1uWDy/AtGkkn7xcgqalcwIQEhFFz5l7UvwUKkqYREREJOnkq2LpCI6lwzebx2A+NI9v2pelkk8WAG7dD6f7jN3ceZByRwNXwiQiIiJJq3gzaPjVo+Xlb+F8cTNTupWnUA7LFCrnbj2g16w9hEZE2SjIJ1PCJCIiIkmv6htQuZ/l39GR8Gs3MgWdZuarFa2jge+/eI+3FhwgKgUObKmESURERJJHw6+gaFPLv8MCYV478tgFMKNHRVwd7QBYc+w6a4+nvCfnlDCJiIhI8ng4sGWuFyzLgVfgl46Uym7PpC7lcbQz81HT4jQq5WXbOGNhMlLb2OQ2FBgYiIeHBwEBAbi7u9s6HBERkdTp/g2YWhfuXbQsF28G7WZzJTCM3JlcEv1wifH7rRYmERERSV4Zc8Arv4Kjm2X5xHJY/2mSJEuJRQmTiIiIJL8cxaHdTDD9m4ps+w4OzLVlRE+khElERERso3A9aDzm0fLyQXB+q83CeRIlTCIiImI7lXpDpb6Wf0dHwMIucPusbWOKhRImERERsa2GX0Gh+pZ/27tARMqbJkUJk4iIiNiWnT20nQ5lOkLvDeBV2tYRPcbe1gGIiIiI4OwOrX+ydRRxUguTiIiISDyUMImIiIjEQwmTiIiISDyUMImIiIjEQwmTiIiISDyUMImIiIjEQwmTiIiISDyUMImIiIjEQwmTiIiISDyUMImIiIjEQwmTiIiISDyUMImIiIjEQwmTiIiISDzsbR1AamIYBgCBgYE2jkRERESe1sPf7Ye/4wmhhOkZBAUFAZA3b14bRyIiIiLPKigoCA8PjwRtazKeJ91KZ6Kjo7l69Spubm6YTKZE3XdgYCB58+bl0qVLuLu7J+q+UxNdBwtdh0d0LSx0HSx0HR7RtbB4mutgGAZBQUHkypULszlhvZHUwvQMzGYzefLkSdJjuLu7p+s3/kO6Dha6Do/oWljoOljoOjyia2ER33VIaMvSQ+r0LSIiIhIPJUwiIiIi8VDClEI4OTkxfPhwnJycbB2KTek6WOg6PKJrYaHrYKHr8IiuhUVyXQd1+hYRERGJh1qYREREROKhhElEREQkHkqYREREROKhhElEREQkHkqYktHEiRPx8fHB2dmZypUrs3v37ifWX7RoEcWKFcPZ2ZnSpUuzatWqZIo0aYwcOZKKFSvi5uZGjhw5aNmyJadOnXriNjNnzsRkMsV4OTs7J1PESWPEiBGPnVOxYsWeuE1aey885OPj89i1MJlM9O/fP9b6aeX9sHnzZpo1a0auXLkwmUwsWbIkxnrDMPjkk0/ImTMnLi4u1KtXj9OnT8e732f9jrG1J12HiIgIhgwZQunSpXF1dSVXrlx069aNq1evPnGfCfl8pQTxvSd69Ojx2Hk1atQo3v2mpfcEEOv3hclkYuzYsXHuM7HeE0qYksnChQt55513GD58OPv376ds2bI0bNiQGzduxFp/+/btdOrUiddee40DBw7QsmVLWrZsydGjR5M58sTj5+dH//792blzJ2vXriUiIoIGDRrw4MGDJ27n7u7OtWvXrK8LFy4kU8RJp2TJkjHOaevWrXHWTYvvhYf27NkT4zqsXbsWgHbt2sW5TVp4Pzx48ICyZcsyceLEWNePGTOG77//nsmTJ7Nr1y5cXV1p2LAhoaGhce7zWb9jUoInXYfg4GD279/Pxx9/zP79+/njjz84deoUzZs3j3e/z/L5Sinie08ANGrUKMZ5/fLLL0/cZ1p7TwAxzv/atWtMnz4dk8lEmzZtnrjfRHlPGJIsKlWqZPTv39+6HBUVZeTKlcsYOXJkrPXbt29vNG3aNEZZ5cqVjb59+yZpnMnpxo0bBmD4+fnFWWfGjBmGh4dH8gWVDIYPH26ULVv2qeunh/fCQ2+99ZZRsGBBIzo6Otb1afH9ABiLFy+2LkdHRxteXl7G2LFjrWX37t0znJycjF9++SXO/Tzrd0xK89/rEJvdu3cbgHHhwoU46zzr5ysliu1adO/e3WjRosUz7Sc9vCdatGhh1KlT54l1Eus9oRamZBAeHs6+ffuoV6+etcxsNlOvXj127NgR6zY7duyIUR+gYcOGcdZPjQICAgDIkiXLE+vdv38fb29v8ubNS4sWLTh27FhyhJekTp8+Ta5cuShQoACdO3fm4sWLcdZND+8FsHxO5s6dS8+ePZ84uXVafD/8v3PnzuHv7x/j/9zDw4PKlSvH+X+ekO+Y1CggIACTyUSmTJmeWO9ZPl+pyaZNm8iRIwdFixbl9ddf5/bt23HWTQ/vievXr7Ny5Upee+21eOsmxntCCVMyuHXrFlFRUXh6esYo9/T0xN/fP9Zt/P39n6l+ahMdHc2gQYOoXr06pUqVirNe0aJFmT59OkuXLmXu3LlER0dTrVo1Ll++nIzRJq7KlSszc+ZMVq9ezaRJkzh37hwvvvgiQUFBsdZP6++Fh5YsWcK9e/fo0aNHnHXS4vvhvx7+vz7L/3lCvmNSm9DQUIYMGUKnTp2eOMHqs36+UotGjRoxe/Zs1q9fz+jRo/Hz86Nx48ZERUXFWj89vCdmzZqFm5sbrVu3fmK9xHpP2D9PsCIJ1b9/f44ePRrvfeSqVatStWpV63K1atUoXrw4P/30E59//nlSh5kkGjdubP13mTJlqFy5Mt7e3vz6669P9ZdSWjVt2jQaN25Mrly54qyTFt8PEr+IiAjat2+PYRhMmjTpiXXT6uerY8eO1n+XLl2aMmXKULBgQTZt2kTdunVtGJntTJ8+nc6dO8f74EdivSfUwpQMsmXLhp2dHdevX49Rfv36dby8vGLdxsvL65nqpyYDBgxgxYoVbNy4kTx58jzTtg4ODpQrV44zZ84kUXTJL1OmTBQpUiTOc0rL74WHLly4wLp16+jVq9czbZcW3w8P/1+f5f88Id8xqcXDZOnChQusXbv2ia1LsYnv85VaFShQgGzZssV5Xmn5PQGwZcsWTp069czfGZDw94QSpmTg6OhI+fLlWb9+vbUsOjqa9evXx/hr+f9VrVo1Rn2AtWvXxlk/NTAMgwEDBrB48WI2bNhA/vz5n3kfUVFRHDlyhJw5cyZBhLZx//59zp49G+c5pcX3wn/NmDGDHDly0LRp02faLi2+H/Lnz4+Xl1eM//PAwEB27doV5/95Qr5jUoOHydLp06dZt24dWbNmfeZ9xPf5Sq0uX77M7du34zyvtPqeeGjatGmUL1+esmXLPvO2CX5PPHe3cXkqCxYsMJycnIyZM2cax48fN/r06WNkypTJ8Pf3NwzDMLp27WoMHTrUWn/btm2Gvb298fXXXxsnTpwwhg8fbjg4OBhHjhyx1Sk8t9dff93w8PAwNm3aZFy7ds36Cg4Ottb573X49NNPjTVr1hhnz5419u3bZ3Ts2NFwdnY2jh07ZotTSBTvvvuusWnTJuPcuXPGtm3bjHr16hnZsmUzbty4YRhG+ngv/L+oqCgjX758xpAhQx5bl1bfD0FBQcaBAweMAwcOGIDx7bffGgcOHLA+/TVq1CgjU6ZMxtKlS43Dhw8bLVq0MPLnz2+EhIRY91GnTh3jhx9+sC7H9x2TEj3pOoSHhxvNmzc38uTJYxw8eDDGd0ZYWJh1H/+9DvF9vlKqJ12LoKAgY/DgwcaOHTuMc+fOGevWrTNeeOEFo3DhwkZoaKh1H2n9PfFQQECAkSFDBmPSpEmx7iOp3hNKmJLRDz/8YOTLl89wdHQ0KlWqZOzcudO6rlatWkb37t1j1P/111+NIkWKGI6OjkbJkiWNlStXJnPEiQuI9TVjxgxrnf9eh0GDBlmvmaenp9GkSRNj//79yR98IurQoYORM2dOw9HR0cidO7fRoUMH48yZM9b16eG98P/WrFljAMapU6ceW5dW3w8bN26M9bPw8Fyjo6ONjz/+2PD09DScnJyMunXrPnZ9vL29jeHDh8coe9J3TEr0pOtw7ty5OL8zNm7caN3Hf69DfJ+vlOpJ1yI4ONho0KCBkT17dsPBwcHw9vY2evfu/Vjik9bfEw/99NNPhouLi3Hv3r1Y95FU7wmTYRjGM7dniYiIiKQj6sMkIiIiEg8lTCIiIiLxUMIkIiIiEg8lTCIiIiLxUMIkIiIiEg8lTCIiIiLxUMIkIiIiEg8lTCIiIiLxUMIkIknKZDI908vHxweA2rVrYzKZOH/+vE3jfxYzZ86McS4ZM2aMsf78+fOYTCZq1679TPv19fWNsd8RI0YkXtAi8lTsbR2AiKRt3bt3f6xs69atnD17lrJly+Lr6xtjXbZs2ZIpsqTz8LycnZ0TZX/NmzfH19eXM2fOsG3btkTZp4g8GyVMIpKkZs6c+VhZjx49OHv2LC1btoyztWT27NkEBweTO3fupA0wCTzpvBLis88+AyzXUgmTiG0oYRKRFClfvny2DkFExEp9mEQkRYqrD9PDfk6RkZF8/vnnFCpUCBcXF4oXL86MGTOs9TZs2MBLL72Eu7s7mTNnplu3bty+fTvWY0VGRjJp0iSqVq2Ku7s7Li4u+Pr68t133xEZGZkk5xcSEsLQoUPx9vbGycmJQoUKMXr0aDQfukjKpBYmEUmV2rdvb02KChYsiJ+fHz179gTAzc2NTp06UaVKFRo2bMiOHTuYM2cO586dY/PmzZhMJut+QkJCaNq0KRs3biRLlixUqVIFZ2dndu3axdtvv83GjRtZvHgxZnPi/X0ZHh5OgwYNOH78OLVr1+bBgwf4+fkxdOhQgoKC+OKLLxLtWCKSOJQwiUiqc+HCBdzc3Dh9+jTZs2cHYOPGjdSpU4cPP/yQ8PBwlixZQtOmTQEIDAykWrVqbN26lU2bNvHSSy9Z9zV48GA2btxIhw4d+Omnn/Dw8AAgKCiIjh07smzZMqZMmUK/fv0SLf4dO3ZQq1Ytzp07h7u7OwB79+6lSpUqjBs3jqFDhz72hJ2I2JZuyYlIqvTdd99ZkyWAl156iXLlynHt2jUaN25sTZYA3N3d6dOnDwB+fn7W8hs3bvDzzz+TN29eZsyYYU2WwNJKNW3aNBwdHZk0aVKixm42m/npp5+syRJAhQoVaNy4McHBwezduzdRjyciz08Jk4ikOg4ODrGOZVSgQAEAGjRoEOe6a9euWcs2bdpEREQEjRo1wsXF5bFtvLy8KFy4MEeOHCEkJCSRogdvb2+KFi36WHmRIkUei1FEUgYlTCKS6nh5eWFnZ/dY+cPbWLENRfBwXVhYmLXsYYfyn3/+Oc6BNI8dO4ZhGNy5cyfR4s+TJ0+s5W5ubo/FKCIpg/owiUiqE18H7KftoB0dHQ1YRtIuW7bsE+s6OTk9XXBPITE7kItI8lDCJCLp1sOWnho1avDDDz/YOBoRScn0Z46IpFsvvfQSdnZ2rFixgoiICFuHIyIpmBImEUm3cufOTc+ePTl//jydOnXi+vXrj9U5c+YMv//+uw2iE5GURLfkRCRdGz9+POfPn+f3339n9erV+Pr6ki9fPh48eMDx48c5c+YMLVq0oE2bNrYOVURsSAmTiKRrLi4u/Pnnn8ybN49Zs2Zx8OBBdu/eTfbs2fH29qZr16507NjR1mGKiI2ZDE1cJCKSKGbOnMmrr77K8OHDGTFiRKrbv4jETS1MIiKJbMmSJZw/fx5nZ2cmT5783Pv75JNPuHjxImfOnEmE6EQkIZQwiYgkskOHDnHo0CFcXV0TJWFatmwZhw4dSoTIRCShdEtOREREJB4aVkBEREQkHkqYREREROKhhElEREQkHkqYREREROKhhElEREQkHkqYREREROKhhElEREQkHkqYREREROKhhElEREQkHv8DA+b7yyCx97IAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for model in models:\n",
    "    time = solutions[model][\"Time [h]\"].entries\n",
    "    voltage = solutions[model][\"Voltage [V]\"].entries\n",
    "    plt.plot(time, voltage, lw=2, label=model.name)\n",
    "plt.xlabel(\"Time [h]\", fontsize=15)\n",
    "plt.ylabel(\"Voltage [V]\", fontsize=15)\n",
    "plt.legend(fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, using `QuickPlot`, we can compare the values of some variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "955c3caa81334b479a57bb7c58b59a48",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=17.0, step=0.17), Output()), _dom_classes=('…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "solution_values = [solutions[model] for model in models]\n",
    "quick_plot = pybamm.QuickPlot(solution_values)\n",
    "quick_plot.dynamic_plot();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we update the current, setting it to be 20 A, we observe a greater discrepancy between the full model and the reduced-order models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "At t = 0.0087774 and h = 1.94278e-62, the corrector convergence failed repeatedly or with |h| = hmin.\n",
      "At t = 0.00384476 and h = 1.59607e-100, the corrector convergence failed repeatedly or with |h| = hmin.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4e115f0bd34b4eb79ff1203344bb04cd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=3415.761176654454, step=34.15761176654454), …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t_eval = np.linspace(0, 3600, 100)\n",
    "for model in models:\n",
    "    solver = pybamm.CasadiSolver()\n",
    "    solutions[model] = solver.solve(model, t_eval, inputs={\"Current function [A]\": 20})\n",
    "\n",
    "# Plot\n",
    "solution_values = [solutions[model] for model in models]\n",
    "quick_plot = pybamm.QuickPlot(solution_values)\n",
    "quick_plot.dynamic_plot();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The relevant papers for this notebook are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n",
      "[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
      "[3] Valentin Sulzer, S. Jon Chapman, Colin P. Please, David A. Howey, and Charles W. Monroe. Faster Lead-Acid Battery Simulations from Porous-Electrode Theory: Part I. Physical Model. Journal of The Electrochemical Society, 166(12):A2363–A2371, 2019. doi:10.1149/2.0301910jes.\n",
      "[4] Valentin Sulzer, S. Jon Chapman, Colin P. Please, David A. Howey, and Charles W. Monroe. Faster Lead-Acid Battery Simulations from Porous-Electrode Theory: Part II. Asymptotic Analysis. Journal of The Electrochemical Society, 166(12):A2372–A2382, 2019. doi:10.1149/2.0441908jes.\n",
      "[5] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pybamm.print_citations()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "19e5ebaa8d5a3277b4deed2928f02ad0cad6c3ab0b2beced644d557f155bce64"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
